快速傅里叶(FFT)算法设计(含程序设计)

快速傅里叶(FFT)算法设计(含程序设计).ppt傅立叶变换将信号从时域转换为频域，可以进行模拟信号的频率分析离散傅立叶变换(DFT)将信号从频域转换为数字(频)域，可以进行数字信号（模拟信号数字化）的频率分析为了实现DFT在计算机上的快速实现，提出了快速离散傅立叶变换(FFT)如何有傅氏变换-DFT-FFT?欧拉公式：=令,称为旋转因子=上式中，k对应数字域，n对应时域另有推导时需用到的公式：1）,lN为l个周期2）,N-m为加上一个周期3),其中4）周期性对称性可约性周期性推导分析若序列x(n)的长度为N，且满足N=2M,(M为自然数)按n的奇偶性把x(n)分解为两个N/2的子序列：x1(r)=x(2r),r=0,1,…,N/2–1x2(r)=x(2r+1),r=0,1,…,N/2–1则x(n)的DFT为：=，k=0,1,…,N/2-1其中X1(k)和X2(k)均以N/2为周期==，k=0,1,…,N/2–1其中公式称为蝶形运算同理，可推出：，k=0,1,…,N/4-1，k=0,1,…,N/4–1……分到最后，k=0，进行蝶形运算的两个输入就是最初输入序列x(n)的其中两个。蝶形分解图示N=8点FFT运算图示N=16点FFT运算图示蝶形运算规律设序列x(n)已经经过时域抽选(倒序)后，存入数组X中。如果蝶形运算的两个输入相距B个点，用原位计算（即计算结果还放在数组的原来位置），则蝶形运算可表示成如下形式：=其中：p=J*2M-L;J=0,1,…,2L-1-1L=1,2,…,M下标L表示第L级运算，XL(J)则表示第L级运算后数组元素X(J)的值。如果要用实数运算完成上述蝶形运算，可由下面的算法进行：(1)(2)设：下标R表示实部下标I表示虚部X’R(J)代表上一次的实数值==(3)(4)(5)(6)(7)==公式(7)(8)(9)主要用于FFT的软件编程由(1)(6)式推导得出由(4)式推导得出由(2)(6)式推导得出由(4)式得出(9)(8)公式中的J就是流程图中公式的变量k，流程图中：N表示阶数，M表示总级数，L表示当前级数，B表示每个蝶形的两个输入数据的间隔，P表示旋转因子指数可以看出，流程图总共有3个循环外循环：次数为级数L的变换范围中循环：为根据当前L求出各个不同的p，循环次数为p的个数2L-1内循环：为每级中每个p对应的蝶形运算个数，循环次数为2M-L内循环中k值每次变换范围(增量)为2L,这是同一级中每个相同的p对应不同蝶形运算的间隔。看图推导软件编程规则：方法一目的是求出旋转因子指数p的变化规律1.当N=8时，第L级共有2L-1个不同的旋转因子。因为N=2M，所以有L=1,2,…,M，即L的最大值为M当L=1时,p=0;(p称为旋转因子指数)当L=2时,p=0,2;k=2(k为p的增量)当L=3时,p=0,1,2,3;k=1当N=16时当L=1时,p=0;当L=2时,p=0,4;k=4当L=3时,p=0,2,4,6;k=2当L=4时,p=0,1,2,3,4,5,6,7;k=12.(j=0,1,2,3,…)(归纳得出：N/k=2L)(L=1,2,3,…表当前级数)(M表总级数)(j=0,1,2…2L-1-1)=p=j*2M-L(变量为j,L)3.每个p对应每级中的运算个数为：第L级中，每个蝶形的两个输入数据相距B=2L-1点同一旋转因子对应着间隔为2L点的2M-L个蝶形看图推导方法二：L=1L=2L=3L=4=M有1个蝶形块，pi=1有2个蝶形块pi=24个蝶形块pi=48个块pi=8pi定义为p的增量反推==pi=2M-L令p=J*pi=J*2M-L(其中J=0,1,2,3,…)，这样写是为了利于软件的循环编程。此时只要求出每级J的个数JTotal即可=JTotal=2L-1得：p=J*2M-L(J=0,1,2,…,2L-1-1)J的总个数JTotal为2L-1，每一级p的总个数也为2L-1外循环次数为级数L中循环为根据当前L求出各个不同的p，循环次数为p的个数2L-1内循环为每级中每个p对应的蝶形运算个数(记为VTotal)，循环次数为2M-L=VTotal=2M-L每个蝶形的两个输入数据间隔（记为INd）：=INd=2L-1同一级中每个相同的p对应蝶形运算的间隔（记为Vd）：=Vd=2L可以看出，为了利于编程，所有的公式推导都往L上靠输入序列倒序的算法设计顺序倒序十进制数I二进制数二进制数十进制数J0000000010011004201001023011110641000011510110156110011371111117顺序与倒序二进制数对照表倒序规律对于N=2M，M位二进制数最高位的权值为N/2，且从左向右二进制位的权值依次为N/4,N/8,…,2,1。因此，最高位加1相当于十进制运算J+N/2。如果最高位是0(JN/2),则直接由J+N/2得下一个倒序值；如果最高位是1(J≥N/2),则要将最高位变成0(J=J-N/2)，次高位加1(J+N/4)。但次高位加1时同样要判断0、1值，如果为0(JN/4),则直接加1(J=J+N/4)，否则先将次高位变成0(J=J-N/4)，再判断下一位，依次类推，直到完成高位加1，适2(1+1=10b)向右进位的运算。I可以从1变换到(N/2-1)，即后半部不用考虑只需前半部和后半部交换后半部不要再重复交换判读各个高位是否为1如果该高位为1，则先减去N/2或N/4、N/8…再判断下个次高位//输入序列倒序软件程序j=N/2;//第0个数(二进制数都为0)和最后一个第N-1个数(二进制数都为1)不需倒序for(i=1;iN-2;i++){if(ij){temp=dataR[i];dataR[i]=dataR[j];dataR[j]=temp;}k=N/2;while(1){if(jk){j=j+k;break;}else{j=j-k;k=k/2;}}}输入序列倒序的算法设计方法二顺序倒序十进制数I二进制数二进制数十进制数J0000000010011004201001023011110641000011510110156110011371111117顺序与倒序二进制数对照表从表格可以看出，所谓倒序只是把数组下标的------最高位与最低位互换次高位与次低位互换…….方法二软件分析：已一个字节(N=256)的倒序为例A[0],A[1],…….,A[255](下标从0000_0000变化到1111_1111)/*定义两个标志位F0，F1*/for(i=0;i≤(255/2);i++)//除2是因为只要判断前半部{j=0;F0=i&0x01;//取最低位F1=i&0x80;//取最高位if(F0)j=j|0x80;//最低位换到最高位if(F1)j=j|0x01;//最高位换到最低位F0=i&0x02;//取次低位F1=i&0x40;//取次高位if(F0)j=j|0x40;//最次位换到次高位if(F1)j=j|0x02;//最次位换到次低位F0=i&0x04;F1=i&0x20;if(F0)j=j|0x20;if(F1)j=j|0x04;F0=i&0x08;F1=i&0x10;if(F0)j=j|0x10;if(F1)j=j|0x08;if(ij)//前半部与后半部交换，相等时无需交换A[i]A[j];}只需单层循环即可，共需要循环（128-2）次算法改进一：前面的算法可以进一步优化为：for(i=0;i≤(255/2);i++){j=0;for(k=0;k4;k++){F0=i&(0x01k);//取最低位F1=i&(0x80k);//取最高位if(F0)j=j|(0x80k);//最低位换到最高位if(F1)j=j|(0x01k);//最高位换到最低位}if(ij)//前半部与后半部交换，相等时无需交换A[i]A[j];}这个算法只是针对8位的，如果是任意N位的应该如何做？这里的N必须满足N=2M算法改进二：针对任意N=2M的情况for(i=0;i(N/2);i++)//或者for(i=0;i≤((N-1)/2);i++){j=0;for(k=0;k(M/2);k++)//当N=256时，M=8{m=0x01k;n=0x80k;F0=i&m;F1=i&n;if(F0)j=j|n;if(F1)j=j|m;}if(ij)//前半部与后半部交换，相等时无需交换A[i]A[j];}FFT软件示例#includemath.h#definePI3.1415926#defineN128#defineM7#defineA0255.0//定义输入波形的幅值voidmain(void){inti,j,k;intp,J,L,B;floatSinIn[N];floatdataR[N],dataI[N];floatdataA[N];floatTr,Ti,temp;//输入波形for(i=0;iN;i++){SinIn[i]=A0*(sin(2*PI*i/25)+sin(2*PI*i*0.4));dataR[i]=SinIn[i];//输入波形的实数部分dataI[i]=0;//输入波形的虚数部分dataA[i]=0;//输出波形的幅度谱为0}//输入序列倒序j=N/2;//第0个数(二进制数都为0)和最后一个第N-1个数(二进制数都为1)不需倒序for(i=1;iN-2;i++){if(ij){temp=dataR[i];dataR[i]=dataR[j];dataR[j]=temp;//因为波形虚数部分都为0,所以不用交换//temp=dataI[i];//dataI[i]=dataI[j];//dataI[j]=temp;}k=N/2;while(1){if(jk){j=j+k;break;}else{j=j-k;k=k/2;}}}//进行FFT//Xr[J]=Xr'(J)+Tr//Xi[J]=Xi'(J)+Ti//Xr[J+B]=Xr'(J)-Tr//Xi[J+B]=Xi'(J)-Ti//(其中Xr'为上一级的Xr,Xi'为上一级的Xi)//其中Tr=Xr'(J+B)cos(2.0*PI*p/N)+Xi'(J+B)sin(2.0*PI*p/N)//Ti=Xi'(J+B)cos(2.0*PI*p/N)-Xr'(J+B)sin(2.0*PI*p/N)for(L=1;L=M;L++)//FFT蝶形级数L从1--M{//计算每个蝶形的两个输入数据相距B=2^(L-1);B=1;i=L-1;while(i){B=B*2;i--;}//或采用运算,即B=2^(L-1);B=(int)pow(2,L-1);//第L级蝶形有pow(2,L-1),即2的L-1次方个蝶形运算和pow(2,L-1)个旋转因子pfor(J=0;J=B-1;J++)//J=0,1,2,...,2^(L-1)-1{//计算p=J*2^(M-L)p=1;i=M-L;while(i){p=p*2;i--;}p=J*p;//第L级蝶形中同一个旋转因子包含多少个蝶形运算//每个蝶形的两个输入数据相距B=2^(L-1)个点//同一旋转因子对应着间隔为2^L点的2^(M-L)个蝶形(2^L=2*B)for(k=J;k=N-1;k=k+2*B){Tr=dataR[k];Ti=dataI[k];temp=dataR[k+B];dataR[k]=dataR[k]+dataR[k+B]*cos(2.0*PI*p/N)+dataI[k+B]*sin(2.0*PI*p/N);dataI[k]=dataI[k]+dataI[k+B]*cos(2.0*PI*p/N)-dataR[k+B]*sin(2.0*PI*p/N);dataR[k+B]=Tr-da