快速傅里叶（FFT）

来源：互联网发布：东北师大网络教育编辑：程序博客网时间：2024/06/06 02:48

快速傅里叶：

更加形象的理解傅里叶变换：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对比结果：

0 0