《信号处理FFT算法.docVIP

实验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=-W FFT把N点DFT运算分解为两组N/2点的DFT运算，然后求和： 其中， 在计算X1(k)与X2(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，蝶形运算可表为： Xm(k)=Xm-1(k)+Xm-1(k+2m-1) WNr Xm(k+2m-1)= Xm-1(k)-Xm-1(k+2m-1) WNr 由于N为已知，所以将r的值确定即可确定WNr。为此，令k=(n2n1n0)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个单元； 存放系数WNr， r=0,1, ,（N/2）-1；需N/2个存储单元； 共计(N+N/2)个存储单元。 5．程序框图 三、实验内容 1．编写一个实现FFT算法的函数。 void fft(COMPLEX A[ ] , int m) 其中，COMPLEX是复数类型的名称，A是要变换的数据，m是FFT点数N的以2为底的幂次(即N=2m )。 编写测试用的主程序对fft函数进行测试，并给出测试结果。 实验代码 #include stdio.h #include stdlib.h #include math.h #define N 1000 #define PI atan(1)*4 typedef 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;isize;i++) { printf(A[%d]:,i);