FFT by C++

単純にC++でFFTを実現し、そこからFFTの速さをより深く理解したいという原理です.しかし初歩的な実現によりFFT,IFFTが実現した.しかし、効率は本当に感動的で、MATLAB 1024*1024サイズのデータFFT+IFFTは0.06秒ぐらいで、私の1024*16サイズのデータは28秒かかります.本当に自分が偽のコードを書くことができないことを感じていますが、私はVectorを使って、コードの実行の効率を全く考慮していません.
兄弟子が書いたFFTを見て、自分は今本当に理解できません.今日はとても軽率な一歩にしましょう.
参考資料:http://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/
ここで、サブデータのサイズと現在のデータサイズの変換に注意します.

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

using namespace std;


const float pi = acos(-1.0);

struct complex
{
    float a, b;
    complex(float __a=0., float __b = 0.): a(__a), b(__b){}
    complex operator + (complex x)
    {
        return complex(a + x.a, b + x.b);
    }
    complex operator * (complex x)
    {
        return complex(a*x.a - b*x.b, a*x.b + b*x.a);
    }
    complex operator / (int n)
    {
        return complex(a / n, b / n);
    }
    complex& operator /= (int n)
    {
        this->a /= n;
        this->b /= n;
        return *this;
    }
};

const int N = 1024*128;

vector fft(vector x)
{
    int len = x.size();
    if (len - (len&(-len))) throw ("fft of recive 2^N!");
    vector ans, even, odd, angle;
    if (len == 1) return x;
    for (int i = 0; i < len; i++)
    {
        if (i & 1) odd.push_back(x[i]);
        else even.push_back(x[i]);
    }

    even = fft(even);
    odd = fft(odd);
    
    for (int i = 0; i < len; i++)
    {
        float a = -2.0*pi*(float)i / (float)len;
        angle.push_back(complex(cos(a), sin(a)));
    }

    for (int i = 0; i < len ; i++)
    {
        x[i] = even[i % (len/2)] + odd[i%(len/2)] * angle[i];
    }

    return x;
}

vector __ifft(vector x)
{
    int len = x.size();
    if (len - (len&(-len))) throw ("fft of recive 2^N!");
    vector ans, even, odd, angle;
    if (len == 1) return x;
    for (int i = 0; i < len; i++)
    {
        if (i & 1) odd.push_back(x[i]);
        else even.push_back(x[i]);
    }

    even = __ifft(even);
    odd = __ifft(odd);

    for (int i = 0; i < len; i++)
    {
        float a = 2.0*pi*(float)i / (float)len;
        angle.push_back(complex(cos(a), sin(a)));
    }

    for (int i = 0; i < len; i++)
    {
        x[i] = (even[i % (len / 2)] + odd[i % (len / 2)] * angle[i]) / (1); //         len/2
    }

    return x;
}

vector ifft(vector x)
{
    x = __ifft(x);
    int N = x.size();
    for (int i = 0; i < N; i++)
    {
        x[i] = x[i] / N;
    }
    return x;
}

int main()
{
    clock_t t = clock();
    vector x;
    for (int i = 0; i < N; i++) x.push_back(complex((float)i, 0.));

    try
    {
        x = fft(x);
    }
    catch (char *s)
    {
        cout << s << endl;
    }


   // x = ifft(x);
    
    
    t = clock();
    cout << t << ' ' << t / CLOCKS_PER_SEC << endl;
    

    //for (int i = 0; i < N; i++)
    //{
    //    cout << i << ": " << x[i].a <0. ? " + ":" ")<< x[i].b<