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 <iostream>
#include <cmath>
#include <vector>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <string>
#include <vector>
#include <ctime>

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<complex> fft(vector<complex> x)
{
    int len = x.size();
    if (len - (len&(-len))) throw ("fft of recive 2^N!");
    vector<complex> 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<complex> __ifft(vector<complex> x)
{
    int len = x.size();
    if (len - (len&(-len))) throw ("fft of recive 2^N!");
    vector<complex> 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<complex> ifft(vector<complex> 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<complex> 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 <<(x[i].b>0. ? " + ":" ")<< x[i].b<<"j" << endl;
    //}

}