4
\$\begingroup\$

I cooked up a FFT code in c++ for calculating with arbitrary data size. The czt function is copied from GNU octave's czt.m. Here's the code,

#define _USE_MATH_DEFINES
#include <vector>
#include <algorithm>
#include <exception>
#include <cmath>
#include <complex>
#include <future>
#include <chrono>
#include <random>
#include <iostream>

std::vector<std::complex<double>> FFT(const std::vector<double>&);
std::vector<std::complex<double>> FFT(const std::vector<std::complex<double>>&);
std::vector<std::complex<double>> CZT(const std::vector<std::complex<double>>&);
std::vector<std::complex<double>> iFFT(const std::vector<std::complex<double>>&);

using namespace std;
using dcvec = vector<complex<double>>;

vector<complex<double>> FFT(const vector<complex<double>>& X)
{
    using namespace std::complex_literals;
    size_t N = X.size();

    if (N == 0) { throw invalid_argument("Empty data"); }

    if (N == 1) { return X; }

    else
    {
        if ((N % 2) != 0) { return CZT(X); }
        else
        {
            vector<complex<double>> y1(N/2), y2(N/2);
            for (size_t n = 0; n < (N / 2); n++)
            {
                y1[n] = X[2 * n];
                y2[n] = X[2 * n + 1];
            }
            y1 = FFT(y1);
            y2 = FFT(y2);

            vector<complex<double>> W_N(N);
            for (size_t n = 0; n < N; n++) { W_N[n] = double(n); }
            for (auto& n : W_N)
            {
                n = exp(-2.0i * M_PI * n / double(N));
            }

            vector<complex<double>> term1(N / 2), term2(N / 2);
            for (size_t n = 0; n < N / 2; n++)
            {
                term1[n] = y1[n] + W_N[n] * y2[n];
                term2[n] = y1[n] + W_N[(N / 2) + n] * y2[n];
            }
            term1.insert(term1.end(), term2.begin(), term2.end());

            return term1;
        }
    }
}

//Same process as GNU Octave's "czt.m"
vector<complex<double>> CZT(const vector<complex<double>>& X)
{
    using namespace std::complex_literals;
    size_t n = X.size();
    complex<double> w = exp(-2.0i * M_PI / double(n));
    double a = 1.0;
    int m = 1 - n;

    vector<complex<double>> chirp(2 * n - 1);
    for (size_t k = 0; k < (2 * n - 1); k++, m++) { chirp[k] = double(m); }
    for (auto& k : chirp) { k = pow(k, 2) / 2.0; }
    for (auto& k : chirp) { k = pow(w, k); }

    size_t N2 = pow(2, ceil(log2(2 * n - 1)));

    vector<complex<double>> xp(n);
    for (size_t k = 0; k < n; k++) { xp[k] = pow(a, -int(k)); }
    for (size_t k = 0; k < n; k++) { xp[k] = X[k] * xp[k] * chirp[k + n - 1]; }
    xp.resize(N2, complex<double>(0, 0));

    vector<complex<double>> ichirp(2 * n - 1);
    for (size_t k = 0; k < (2 * n - 1); k++) { ichirp[k] = 1.0 / chirp[k]; }
    ichirp.resize(N2, complex<double>(0, 0));

    auto fut_xp = async(launch::deferred, (dcvec(*) (const dcvec&)) &FFT, xp);
    auto fut_ichirp = async(launch::deferred, (dcvec(*) (const dcvec&)) &FFT, ichirp);
    xp = fut_xp.get();
    ichirp = fut_ichirp.get();
    for (size_t k = 0; k < N2; k++) { xp[k] = xp[k] * ichirp[k]; }
    xp = iFFT(xp);

    vector<complex<double>> ret(n);
    for (size_t k = 0; k < n; k++) { ret[k] = xp[k + n - 1] * chirp[k + n - 1]; }

    return ret;
}

vector<complex<double>> FFT(const vector<double>& X)
{
    size_t n = X.size();
    vector<complex<double>> Y(n);
    for (size_t i = 0; i < n; i++)
    {
        Y[i] = complex<double>(X[i], 0);
    }

    return FFT(Y);
}

vector<complex<double>> iFFT(const vector<complex<double>>& X)
{
    size_t N = X.size();
    vector<complex<double>> ret = X;
    for (auto& n : ret) { n = conj(n); }
    ret = FFT(ret);
    for (auto& n : ret) { n = conj(n) / double(N); }

    return ret;
}

And here's the driver code,

int main()
{
    size_t n = 910;
    vector<double> A(n, 0);
    mt19937 gen(10);
    uniform_real_distribution<double> rng(-10.0, 10.0);

    for (size_t i = 0; i < n; i++)
    {
        A[i] = rng(gen);
    }
    
    auto start = chrono::high_resolution_clock::now();
    auto B = FFT(A);
    auto end = chrono::high_resolution_clock::now();
    cout << chrono::duration<double>(end - start) << endl;

    return 0;
}

You can copy paste everything into one file and it should work. The results are actually correct and it can work with any size of the data. Just play around with the numbers for size_t n = #### in main function.

But the performance is not really as good as I wanted it to be. Now, on my computer, it takes roughly 60~80ms for a data size of 910. The same data can be fft-ed by scipy under 2ms. So, I am guessing my code can be optimized to achieve at least <10ms timing.

Please, give me advice on how to improve the timing. Any advice is appreciated. Thank you for your time.

Note: I would prefer not to use third-party libraries due to circumstances. This is not a classwork/HW/assignment.

\$\endgroup\$
4
  • 1
    \$\begingroup\$ Request for clarification. What is the point of for (auto& n : W_N){ n = exp(-2.0i * M_PI * n / double(N)); }? What is the point of for (auto& k : chirp) { k = pow(k, 2) / 2.0; } and for (auto& k : chirp) { k = pow(w, k); }? \$\endgroup\$
    – vnp
    Jun 3, 2022 at 21:41
  • \$\begingroup\$ The first one is the complex coefficients in Cooley-Tukey algorithm. The second one (chirp), I don't know. czt.m had it, I have it. I still don't understand these math for chirp-z transform. \$\endgroup\$ Jun 4, 2022 at 0:02
  • \$\begingroup\$ Are you compiling with optimizations? The code is already <10ms as is. You don't have to store W_N in an array, just compute these values as they're needed. You can also notice that the second half of W_N is equal to the first half negated and halve the number of exponentiations. Then you can observe that each element of W_N differs from the previous one by a constant factor and replace most exponentiations with multiplications. Then you can replace the only remaining exponentiation with a table lookup, this is enough to make the function take <1ms on my machine. \$\endgroup\$ Jun 4, 2022 at 4:38
  • \$\begingroup\$ @thedefault. compiling with optimization didn't help much. Actually, given a data size of 1024, it finishes in ~10ms. But for some reason the two async calls at the end of czt takes ~20ms, like they are executing sequentially. As for the rest, I am trying to reduce the unnecessary containers like W_N. I don't really understand the math of your comment's last part but I will try. \$\endgroup\$ Jun 4, 2022 at 14:33

1 Answer 1

1
\$\begingroup\$

As it turned out, the issue was the compiler. I was using msvc which does bound check for vector subscript. Since the code has a lot of computation on vector using subscript, the bound check consumes a lot of time as well. I assume using iterator instead would speed it up. However, compiling in g++ with -O3 flag gives a much better performance (~6ms). I did make some more adjustments to the code as well, thanks to @the default for his tips. I will put the code here for anyone interested.

#define _USE_MATH_DEFINES
#include <vector>
#include <algorithm>
#include <exception>
#include <cmath>
#include <complex>
#include <future>
#include <iostream>
#include <chrono>
#include <random>
#include <fstream>

using namespace std;
using dcvec = vector<complex<double>>;

vector<complex<double>> FFT(const vector<double>&);
vector<complex<double>> FFT(const vector<complex<double>>&);
vector<complex<double>> CZT(const vector<complex<double>>&);
vector<complex<double>> iFFT(const vector<complex<double>>&);

vector<complex<double>> FFT(const vector<complex<double>>& X)
{
    using namespace std::complex_literals;
    size_t N = X.size();

    if (N == 0) { throw invalid_argument("Empty data"); }

    if (N == 1) { return X; }

    else
    {
        if ((N % 2) != 0)
        {
            return CZT(X);
        }
        else
        {
            vector<complex<double>> y1(N / 2), y2(N / 2);
            for (size_t n = 0; n < (N / 2); n++)
            {
                y1[n] = X[2 * n];
                y2[n] = X[2 * n + 1];
            }
            y1 = FFT(y1);
            y2 = FFT(y2);

            vector<complex<double>> term1(N / 2), term2(N / 2);
            complex<double> W_N = 1.0;
            complex<double> fact = exp(-2.0i * M_PI / double(N));

            term1[0] = y1[0] + y2[0];
            term2[0] = y1[0] - y2[0];

            for (size_t n = 1; n < N / 2; n++)
            {
                W_N = W_N * fact;
                term1[n] = y1[n] + W_N * y2[n];
                term2[n] = y1[n] + (-W_N) * y2[n];
            }
            term1.insert(term1.end(), term2.begin(), term2.end());

            return term1;
        }
    }
}

//Same process as GNU Octave's "czt.m"
vector<complex<double>> CZT(const vector<complex<double>>& X)
{
    using namespace std::complex_literals;
    size_t n = X.size();
    size_t N = 2 * n - 1;
    complex<double> w = exp(-2.0i * M_PI / double(n));
    double a = 1.0;
    int m = 1 - n;

    vector<complex<double>> chirp(N);
    for (size_t k = 0; k < N; k++, m++) { chirp[k] = pow(w, (pow(double(m), 2) / 2.0)); }

    size_t N2 = pow(2, ceil(log2(N)));

    vector<complex<double>> xp(n);
    for (size_t k = 0; k < n; k++) { xp[k] = X[k] * pow(a, -int(k)) * chirp[k + n - 1]; }
    xp.resize(N2, complex<double>(0, 0));

    vector<complex<double>> ichirp(N);
    transform(chirp.begin(), chirp.end(), ichirp.begin(), [](const complex<double>& k)->complex<double> {return (1.0 / k); });
    ichirp.resize(N2, complex<double>(0, 0));

    auto fut_xp = async(launch::deferred, (dcvec(*) (const dcvec&)) & FFT, xp);
    auto fut_ichirp = async(launch::deferred, (dcvec(*) (const dcvec&)) & FFT, ichirp);
    xp = fut_xp.get();
    ichirp = fut_ichirp.get();
    for (size_t k = 0; k < N2; k++) { xp[k] = xp[k] * ichirp[k]; }
    xp = iFFT(xp);

    vector<complex<double>> ret(n);
    for (size_t k = 0; k < n; k++) { ret[k] = xp[k + n - 1] * chirp[k + n - 1]; }

    return ret;
}

vector<complex<double>> FFT(const vector<double>& X)
{
    size_t n = X.size();
    vector<complex<double>> Y(n);
    for (size_t i = 0; i < n; i++)
    {
        Y[i] = complex<double>(X[i], 0);
    }

    return FFT(Y);
}

vector<complex<double>> iFFT(const vector<complex<double>>& X)
{
    size_t N = X.size();
    vector<complex<double>> ret = X;
    for (auto& n : ret) { n = conj(n); }
    ret = FFT(ret);
    for (auto& n : ret) { n = conj(n) / double(N); }

    return ret;
}

int main()
{
    size_t n = 910;
    vector<double> A(n, 0);
    mt19937 gen(10);
    uniform_real_distribution<double> rng(-10.0, 10.0);

    for (size_t i = 0; i < n; i++)
    {
        A[i] = rng(gen);
    }

    auto start = chrono::high_resolution_clock::now();
    auto B = FFT(A);
    auto end = chrono::high_resolution_clock::now();
    cout << chrono::duration<double>(end - start).count() << endl;

    ofstream FILE("Data.txt");

    for (const auto &n : A) { FILE << n << "\t"; }

    FILE.close();

    FILE.open("Result.txt");

    for (const auto &n : B) {
        FILE << n.real();
        if (n.imag() < 0) { FILE << n.imag() << "j" << "\t"; }
        else { FILE << "+" << n.imag() << "j" << "\t"; }
    }

    FILE.close();

    return 0;
}

The file output is to compare the result with scipy using this script,

import time
from scipy.fft import fft
import numpy as np

file = open("Data.txt", 'r')
line = file.readline()
file.close()
line = [float(x) for x in line.split("\t") if x!=""]
x = np.asarray(line)
t1 = time.time_ns()
print(t1)
y = fft(x)
t2 = time.time_ns()
print(t2)
t = t2-t1
print(t)

file = open("Result.txt")
line = file.readline()
line = [complex(x) for x in line.split("\t") if x!=""]
y2 = np.asarray(line)

t = abs(y-y2)
for i in t:
    if i > 0.001:
        print("{:1.0}".format(i))

If anyone has further improvement to propose, feel free to add an answer.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.