FFT算法设计(含程序设计).ppt

《FFT算法设计(含程序设计).ppt》由会员分享，可在线阅读，更多相关《FFT算法设计(含程序设计).ppt（37页珍藏版）》请在课桌文档上搜索。

1、第6讲 FFT算法设计,傅立叶变换将信号从时域转换为频域，可以进行模拟信号的频率分析离散傅立叶变换(DFT)将信号从频域转换为数字(频)域，可以进行数字信号（模拟信号数字化）的频率分析为了实现DFT在计算机上的快速实现，提出了快速离散傅立叶变换(FFT),如何有傅氏变换-DFT-FFT?,欧拉公式：=令,称为旋转因子=上式中，k对应数字域，n对应时域,另有推导时需用到的公式：1）,l N为l个周期 2）,N-m为加上一个周期3),其中4）,周期性,对称性,可约性,周期性,推导分析,若序列x(n)的长度为N，且满足N=2M,(M为自然数)按n的奇偶性把x(n)分解为两个N/2的子序列：x1(r)

2、=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个点，用原位计算（即计算结果还放在数组的原来位置），则蝶形运算可表

3、示成如下形式：=其中：p=J*2M-L;J=0,1,2L-1-1 L=1,2,M 下标L表示第L级运算，XL(J)则表示第L级运算后数组元素X(J)的值。如果要用实数运算完成上述蝶形运算，可由下面的算法进行：,(1),(2),设：下标R表示实部 下标I表示虚部 XR(J)代表上一次的实数值=,(3),(4),(5),(6),(7),=公式(7)(8)(9)主要用于FFT的软件编程,由(1)(6)式推导得出,由(4)式推导得出,由(2)(6)式推导得出,由(4)式得出,(9),(8),公式中的J就是流程图中公式的变量k，流程图中：N表示阶数，M表示总级数，L表示当前级数，B表示每个蝶形的两个输入

4、数据的间隔，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

5、;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=1,2.(j=0,1,2,3,)(归纳得出：N/k=2L)(L=1,2,3,表当前级数)(M表总级数)(j=0,1,22L-1-1)=p=j*2M-L(变量为j,L),3.每个p对应每级中的运算个数为：第L级中，每个蝶形的两个输入数据相距B=2L-1点 同一旋转因子对应着间隔为2L点的2M-L个蝶形,看图推导方法二：,L=1,L=2,L=3,L=4=M,有1个蝶形块，pi=1,有2个蝶形块pi=2,4个蝶形块pi=4,8个块pi

6、=8,pi定义为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

7、对应蝶形运算的间隔（记为Vd）：,=Vd=2L,可以看出，为了利于编程，所有的公式推导都往L上靠,输入序列倒序的算法设计,顺序与倒序二进制数对照表,倒序规律,对于N=2M，M位二进制数最高位的权值为N/2，且从左向右二进制位的权值依次为N/4,N/8,2,1。,因此，最高位加1相当于十进制运算J+N/2。如果最高位是0(JN/2),则直接由J+N/2得下一个倒序值；,如果最高位是1(JN/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)，再判断下一位

8、，依次类推，直到完成高位加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;i N-2;i+)if(i j)temp=dataRi;dataRi=dataRj;dataRj=temp;k=N/2;while(1)if(j k)j=j+k;break;else j=j-k;k=k/2;,输入序列倒序

9、的算法设计方法二,顺序与倒序二进制数对照表,从表格可以看出，所谓倒序只是把数组下标的-最高位与最低位互换次高位与次低位互换.,方法二软件分析：已一个字节(N=256)的倒序为例,A0,A1,.,A255(下标从0000_0000变化到1111_1111)/*定义两个标志位F0，F1*/for(i=0;i(255/2);i+)/除2是因为只要判断前半部 j=0;F0=i/最次位换到次低位,F0=i,只需单层循环即可，共需要循环（128-2）次,算法改进一：,前面的算法可以进一步优化为：for(i=0;i(255/2);i+)j=0;for(k=0;kk);/取最高位 if(F0)j=j|(0 x

10、80k);/最低位换到最高位 if(F1)j=j|(0 x01k);/最高位换到最低位 if(ij)/前半部与后半部交换，相等时无需交换 Ai Aj;,这个算法只是针对8位的，如果是任意N位的应该如何做？这里的N必须满足N=2M,算法改进二：针对任意N=2M的情况,for(i=0;ik;F0=i,FFT软件示例,#include#define PI 3.1415926#define N 128#define M 7#define A0 255.0/定义输入波形的幅值void main(void)int i,j,k;int p,J,L,B;float SinInN;float dataRN,da

11、taIN;float dataAN;float Tr,Ti,temp;/输入波形for(i=0;i N;i+)SinIni=A0*(sin(2*PI*i/25)+sin(2*PI*i*0.4);dataRi=SinIni;/输入波形的实数部分dataIi=0;/输入波形的虚数部分dataAi=0;/输出波形的幅度谱为0,/输入序列倒序j=N/2;/第0个数(二进制数都为0)和最后一个第N-1个数(二进制数都为1)不需倒序for(i=1;i N-2;i+)if(i j)temp=dataRi;dataRi=dataRj;dataRj=temp;/因为波形虚数部分都为0,所以不用交换/temp=d

12、ataIi;/dataIi=dataIj;/dataIj=temp;k=N/2;while(1)if(j k)j=j+k;break;else j=j-k;k=k/2;,/进行FFT/XrJ=Xr(J)+Tr/XiJ=Xi(J)+Ti/XrJ+B=Xr(J)-Tr/XiJ+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

13、/计算每个蝶形的两个输入数据相距 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)个点/同一旋转因子对应着间隔为2L点的2(M-L)个蝶形(2L=2*B)

14、for(k=J;k=N-1;k=k+2*B)Tr=dataRk;Ti=dataIk;temp=dataRk+B;dataRk=dataRk+dataRk+B*cos(2.0*PI*p/N)+dataIk+B*sin(2.0*PI*p/N);dataIk=dataIk+dataIk+B*cos(2.0*PI*p/N)-dataRk+B*sin(2.0*PI*p/N);dataRk+B=Tr-dataRk+B*cos(2.0*PI*p/N)-dataIk+B*sin(2.0*PI*p/N);dataIk+B=Ti-dataIk+B*cos(2.0*PI*p/N)+temp*sin(2.0*PI*p

15、/N);,/求出幅度频率谱/因为从0-N是0-2PI范围,所以只要求出0-N/2即可/幅频对应的位置可由128*(1/2PI)=(128/25)*(1/x)求出/-/|/(采样频率)1屏总个数 周期(输入频率)1屏总个数 输入频率周期/得出输入频率x=2PI/25/对应的幅频值的波形位置=128/25=5.12(因为2PI对应点的位置为128,PI对应点的位置为64)for(i=0;iN/2;i+)dataAi=sqrt(dataRi*dataRi+dataIi*dataIi);/如果要算相位，则i=arctan(dataIi/dataRi)/频谱的平方称为功率谱，记为：Pi=dataRi*dataRi+dataIi*dataIiwhile(1);/break point,