Ways to optimize my FFT impementation's performance

Question

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.

Answer 1

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.