信号处理 FFT算法

实验2 基2时域抽选的FFT 程序设计与调试一、实验目的掌握信号处理，尤其是数字信号处理的基本原理和方法。

要求能通过实验熟练掌握基2时域抽选的快速傅立叶变换算法（FFT ）的基本原理，了解二维及多维快速傅立叶变换算法。

二、实验原理1．复数类型对于FFT 算法涉及的复数运算，使用自定义的COMPLEX 来定义复数类型，其使用方法与常规类型（如int,float,double ）相似。

typedef struct { float real, imag; } COMPLEX; 2．FFT 基本原理FFT 改进了DFT 的算法，减少了运算量，主要是利用了旋转因子W 的两个性质：（a ）W 的周期性：W = W （b) W 的对称性：W =-WFFT 把N 点DFT 运算分解为两组N/2点的DFT 运算，然后求和：)()()(21k X W k X k X k N +=1,,1,0 ),()()2(221-=-=+N kN k k X W k X N k X 其中，∑∑∑∑-=-=-=-=+====1122111122222222)12()()()2()()(N N NN N N N N r rkr rk r rkr rkW r x W r xk X Wr x Wr x k X在计算X 1(k)与X 2(k)时，仍利用上述公式，把它们看成是新的X(k)。

如此递归下去，便是FFT 算法。

3．蝶形运算从基2时域抽选FFT 运算流图可知：① 蝶形两节点的距离为2m-1，其中，m 表示第m 列，且m =1,… ,L 。

例如N=8=23， 第一级(列)距离为21-1=1， 第二级(列)距离为22-1=2， 第三级(列)距离为23-1=4。

② 考虑蝶形运算两节点的距离为2m-1，蝶形运算可表为： X m (k)=X m-1(k)+X m-1(k+2m-1) W N rX m (k+2m-1)= X m-1(k)-X m-1(k+2m-1) W N r由于N 为已知，所以将r 的值确定即可确定W N r 。

为此，令k=(n 2n 1n 0)2 ,再将k 左移(L-m)位，右边位置补零，就可得到(r)2 的值，即(r)2 =(k)22L-m 。

例如 N=8=23(1) k=2 , m=3 的r 值,∵k=2=(010)2 左移L-m=3-3=0 ,∴ r=(010)2=2; (2) k=3 ,m=3 的r 值;k= 3=(011)2,左移0位,r=3；(3) k=5 ,m=2的值;k=5=(101)2 左移L-m=1位 r=(010)2 =2。

4．存储单元存输入序列 x(n)，n=0,1, ,N-1，计N 个单元；存放系数W N r ， r=0,1, ,（N/2）-1；需N/2个存储单元； 共计(N+N/2)个存储单元。

5．程序框图开 始送入x (n )，MN ＝2 M倒 序L ＝1 , MＪ＝0 , B £ 1P ＝2 M £L Jk ＝ J , N £1 , 2LpNpNW B k X k X B k X W B k X k X k X )()()()()()(+-⇐+++⇐输出结 束B 2 L £1三、实验内容1．编写一个实现FFT 算法的函数。

void fft(COMPLEX A[ ] , int m) 其中，COMPLEX 是复数类型的名称，A 是要变换的数据，m 是FFT 点数N 的以2为底的幂次(即N=2m )。

2．编写测试用的主程序对fft 函数进行测试，并给出测试结果。

四、实验代码#include <stdio.h> #include <stdlib.h>#include <math.h>#define N 1000#define PI atan(1)*4typedef struct{double real;double imag;}complex;void fft(complex A[],int m); //快速傅里叶变换void initW();void change();void add(complex,complex,complex *); //复数加法void mul(complex,complex,complex *); //复数乘法void sub(complex,complex,complex *); //复数减法void divi(complex,complex,complex *); //复数除法void output();complex A[N], *W; //输出序列的值complex up,down,product;int size=0; //序列的长度int m;void main(){fft(A,m);output();}//进行基-2 FFT运算void fft(complex x[],int m){ int i=0,j=0,k=0,l=0;printf("Please input m:");scanf("%d",&m);size=pow(2,m);printf("Please input %d datas in A[%d]:\n",size,size);for(i=0;i<size;i++){printf("A[%d]:",i);scanf("%lf %lf",&A[i].real,&A[i].imag);}initW();change();for(i=0;i<log(size)/log(2);i++){l=1<<i;for(j=0;j<size;j+=2*l){for(k=0;k<l;k++){mul(A[j+k+l],W[size*k/2/l],&product);add(A[j+k],product,&up);sub(A[j+k],product,&down);A[j+k]=up;A[j+k+l]=down;}}}}void initW() //初始化，得到W算子的值{int i;W=(complex *)malloc(sizeof(complex)* size);for(i=0;i<size;i++){W[i].real=cos(2*PI/size*i);W[i].imag=-1*sin(2*PI/size*i);}}void change() //对时间序列进行倒序{complex temp;unsigned int i=0,j=0,k=0;double t;for(i=0;i<size;i++){k=i;j=0;t=(log(size)/log(2));while((t--)>0){j=j<<1; //左移一位j|=(k&1); //与1进行与运算k=k>>1; //右移一位}if(j>i){temp=A[i];A[i]=A[j];A[j]=temp;}}}void output(){int i;printf("The results are:\n");for(i=0;i<size;i++){printf("%.4f",A[i].real);if(A[i].imag>=0.0001)printf("+%.4fj\n",A[i].imag);else if(fabs(A[i].imag)<0.0001)printf("\n");elseprintf("%.4fj\n",A[i].imag);}}void add(complex a,complex b,complex *c){c->real=a.real+b.real;c->imag=a.imag+b.imag;}void mul(complex a,complex b,complex *c){c->real=a.real*b.real - a.imag*b.imag;c->imag=a.real*b.imag + a.imag*b.real;}void sub(complex a,complex b,complex *c){c->real=a.real-b.real;c->imag=a.imag-b.imag;}void divi(complex a,complex b,complex *c){c->real=(a.real*b.real+a.imag*b.imag)/(b.real*b.real+b.imag*b.imag);c->imag=(a.imag*b.real-a.real*b.imag)/(b.real*b.real+b.imag*b.imag); }五、实验截图。