快速傅里叶（FFT）

快速傅里叶：

更加形象的理解傅里叶变换：http://zhuanlan.zhihu.com/wille/19763358

大概了解之后：https://www.zhihu.com/question/19714540?rf=20254033

从傅里叶级数到傅里叶变换：

（照片太大，只能裁剪为两张）

刨根问底的同学：https://www.zhihu.com/question/21665935

雷德算法：

输出序列是按自然顺序排列的，而输入序列的顺序则是“比特反转”方式排列的。也就是说，将序号用二进制表示，然后将二进制数以相反方向排列，再以这个数作为序号。如011变成110，那么第3个输入值和第六个输入值就要交换位置了。

详见http://blog.csdn.net/axiqia/article/details/50830435

思路：

代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#define PI 3.1415926   //定义圆周率值
#define N  16 //变换点数

using namespace std;

typedef struct compx
{
    float real;
    float img;
}compx;


//复数乘法
compx Multi(compx a, compx b)
{
    compx c;
    c.real = a.real*b.real - a.img*b.img;
    c.img = a.real*b.img + a.img*b.real;
    return c;
}
//复数加法
compx Add(compx a, compx b)
{
    compx c;
    c.real = a.real + b.real;
    c.img = a.img + b.img;
    return c;
}
//复数减法
compx Sub(compx a, compx b)
{
    compx c;
    c.real = a.real - b.real;
    c.img = a.img - b.img;
    return c;
}

//快速傅里叶
void FFT(compx data[])
{
    int k, i, j = 0;
    compx u, w, t;
    //自然顺序变成倒位序

    for(i = 0; i < N-1; i++)
    {
        if(i < j)
        {
            t.real = data[i].real;
            t.img = data[i].img;
            data[i].real = data[j].real;
            data[i].img = data[j].img;
            data[j].real = t.real;
            data[j].img = t.img;
        }
        k = N/2;

        //实现进位
        while(k <= j)   //左边最高位为1
        {
            j = j-k;    //最高位1变为0
            k /= 2;
        }
        j += k;

    }

    //计算蝶形级数m
    int f = N, m;
    for(m = 1; (f=f/2) != 1; m++);


    //控制级数
    for(int n = 1; n <= m; n++)
    {
        int a = 2<<(n-1);      //相邻蝶形结的距离.2的n次方
        int b = a/2;          //相邻运算点的距离
        //printf("a = %d b = %d\n", a,b);
        u.real = 1.0;       //蝶形运算的系数
        u.img = 0.0;
        w.real = cos(PI/b);//当前系数与上一系数的商
        w.img = -sin(PI/b);

        int x1, x2;     //参加蝶形运算的两个节点
        for(j = 0; j < b; j++)
        {
            //printf("b = %d\t%lf\t%lf\n",b, u.real, u.img);
            for(x1 = j; x1 < N; x1 += a)
            {
                x2 = x1 + b;

                t = Multi(data[x2], u);
                data[x2] = Sub(data[x1], t);
                data[x1] = Add(data[x1], t);
                //printf("x1 = %d, x2 = %d\n", x1, x2);
            }
            u = Multi(u, w);
        }
    }

}
int main()
{
    compx data[N];
    for(int i = 0; i < N; i++)
    {
        data[i].real = i;
        data[i].img = 0;
    }

    FFT(data);
    for(int i = 0; i < N; i++)
        printf("%lf \t%lf\n", data[i].real, data[i].img);

    return 0;
}

运行结果：

MATLAB对比结果：