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

1、第6讲FFT算法设计傅立叶变换将信号从时域转换为频域，可以进行 模拟信号的频率分析离散傅立叶变换(DFT)将信号从频域转换为数字 (频)域，可以进行数字信号(模拟信号数字化) 的频率分析为了实现DFT在计算机上的快速实现，提出了快速离散傅立叶变换(FFT)如何有傅氏变换.DFT- FFT?e j(,) = cosco j sin欧拉公式：尸0)=匚 fe-jwtdt = f(t)edt = f(t)edt-0 knN7.j 二= X(k) = A(/z)e “ 仏=0丄 2,Nfi=0令必=不F称为旋转因子N_X 伙)=工 x(n)Wn, k = 0,1,2,“ 1上式中，k对应数字域，n对应

2、时域另有推导时需用到的公式：1)W+IN = e小,1 N为I个周期jm=e N = W；17、-7*(-w)Z) 乂丁 =e n=WmN2兀3) % 2 一I对称性周期性j聲e，Nm为加上一个周期|周期性僅韦,N、=-4）叽一F-jm一-Je N *e N 2 =，其中 ej7T cos(tt) Jsin() = -12兀-J*. N _一 n可约哲推导分析若序列x(n)的长度为N,且满足N = 2M,(M为自然数) 按n的奇偶性把x(n)分解为两个N/2的子序列： x1(r)=x(2r)z r=0zlz.vN/2 - 1 x2(r)=x(2r+l)z r=Oflz.vN/2 - 1则x(n

3、)的DFT为：X(R)=工.心)歐” + 工，匸偶数，匸奇数gN/2-1N/2-I=x(2r)Wkr + 为班 2厂+ l)W?m)r=()r=()N/2-1N/2-1=石(厂)-勺什)吒N/2-1工 xgWd + W： 22(门咋2N/2-1r=0r=0r=0r=0X(k) = XSk)+WX2(k),k=0,lz.zN/2 1其中XJk）和X2（k）均以N/2为周期NN 梯NX(k +于)=X|伙+亍)+叫2X?伙+刁)x.（0）X3(0)Xi(l)X(2)X.(3)&(0)M2)心3)N/4点DFTN/4点DFTN/4点DFTDFT蝶形分解图示x(0)x(4) x(2) x(6) x(l

4、) x(5) M3) x(7)X4(lW3同理，可推出：x )=X?伙)+ W爲总4仏)n， k=Ozlz.,N/4 - 1X + -)=X.(k)-W,2X4(k)4X?伙） X5伙）十w爲2乂6伙）x十牛=X5(k)-W爲2X9 k=0几,N/4 - 1分到最后，k=0,进行蝶形运算的两个输入就是最初输入序 列x（n）的其中两个。X(0) x(l) X X X(4) X X(6) X(7)N=8点FFT运算图示x(0)4)3M2)%(6)f竺)x(l) M也x(5)朋x(3) 3也x(7)皿辺A(0)A(0)A(0)A(l)A(2)A(2)AAAAA(7)A(l)A(6)妙 x(o) AO

5、l x(l)X 迥X A(4 x X(5) X(6) &Zlx(7)N = 16点FFT运算图示蝶形运算规律设序列x（n）已经经过时域抽选（倒序）后，存入数组X中。如 果蝶形运算的两个输入相距B个点，用原位计算（即计算结果 还放在数组的原来位置），则蝶形运算可表示成如下形式：X伙）=X（伙）十呎X2伙）X（Z普）=X（約一比兀（幻X,（丿）=X-（丿）+ X丿 + B）Wf；X/（y + ）v=X/gz） X, + （并其中：p = J*2ML; J = 0,lf.f2L1-l L=1Z2,.,M下标L表示第L级运算，XL（J）则表示第L级运算后数组元素 XQ）的值。如果要用实数运算完成上述蝶

6、形运算，可由下面的算法进行:设：T = XZ=TR+jTl （3）下标R表示实部X/ 1（J） = X/?（/） + jX,J） （4）下标1 表示虚部 X）Q）代表上一次的实数值X-(八 3) = (八 8) + JX/(J + B) (5) X_G/ + 3)W/ =XQ + 3)fX/G/ + 3)W/= lXK(J+B) + jX/(J + B)|(cos 等 p-./sin 等/?)= XR(J + B)cos p+ X,(J + B)sin p+ jX,(J + B)cos p XN(J + B)sin p NNNN= Tr5(6)Tr =Xr(J + B)cos# p + x；(

7、丿 + B)sin 千 p 7； = X/ (丿 + B)cos# p 一 Xr(丿 + 5)sin pXL(J) = XR(J) + jX/(J)= XL_SJ)+TR+jT由(6)式推导得出=XRU) + jXm+TR +丿刀由式推导得出X.(丿 + 3) = X- (J) 一 (7； + JTf)由(2)(6)式推导得出= X、RQ) + jX,J)-(jR+ j)由(4)式得出(8)(9)= 厂应心乂斥+几X/(J) = X；(J) + T/Xr(J + B) = Xk(J)-TrX/(J + B) = X/(J)-T/送入 x(n), N, M开始输入序列倒序公式(7)(8)(9)主

8、要用于FFT的软件编程Tr = X,八 )cos# “ + X/(丿 + B)sin #T, = X / (丿 + B)cos- p - X N(J + B)sin 二 pXr(J) = Xr(J) + TrX/(J) = X；(J)+T/Xr(J+B) = Xr(J) TrB=2A(L-1)P= J*2A(M-L)X伙)v= X(A)+ X(Ar + B)W X(k + B) pi = 2M_L令p=J*pi=J*2ML （其中J=0,1,2,3,），这样写是为了利于软件的循环编程。此时只要求出每级J的个数J”即可当L=lfl寸，J当厶=2时，J当厶=3时，片罰TotalTotal= 1 =

9、 2= 2 = 2*= 4 = 2?= JTotal=2L1% = 4时，JTolaI =8 = 2得：p = J*2M L(J=0,1,2,.,2-1)J的总个数J“讪为2, 每一级P的总个数也为2外循环次数为级数L中循环为根据当前L求出各个不同的p,循环次数为p的个 数2-丄内循环为每级中每个p对应的蝶形运算个数（记为）， 循环次数为2丄当乙=1时，VTixal = 8 = 23Total当乙=2时，Vt如=4 = 2? 当乙=3时，.=2 = 2*当无=4时，J侖=1 = 2每个蝶形的两个输入数据间隔（记为INd）:当厶=1吋，LV = 1 = 2 当 = 2H寸，IN(I =2 = 2

10、 当厶=3时，Wd =4=22 当厶=4时，Wtl = 8 = 23同一级中每个相同的P对应蝶形运算的间隔（记为VQ ：= Vd = 2L当L = W寸,Vz = 2 = 2, 当L = 2B寸,Vz = 4 = 2?当L = 3H寸，V z = 8 = 23 当厶=4(甘，V =16=24可以看出，为了利于编程，所有的公式推导都往L上靠输入序列倒序的算法设计顺序倒序十进制数I二进制数二进制数十进制数J000000001001100423010011010110264567100 1011110111(0011010111111537x(0) x(l) x(2) x(3) x(4) x(5)

11、x(6) x(7)A(0) A A(2) A(3) A(4) A(5) A(6) A(7)A(0) A(l) A(2) A(3) A(4) A(5) A(6) A(7) x(0) x(4) x(2) x(6) x(l) x(5) x(3) x(7)倒序规律对于N=2M, M位二进制数最高位 的权值为N/2,且从左向右二进制 位的权值依次为N /4,N/8,.,2,1。因此，最高位加1相当于十进制运 算J+N/2o如果最高位是0(JvN/2), 则直接由J+N/2得下一个倒序值； 如果最高位是1(JNN/2),则要将最高位变成O(Jv二JN/2),次高位 加1(J+N/4)。但次高位加1时同样要

12、判断0、1值，如果为0(JN/4), 则直接加1(J=J+N/4),否则先将次高位变成0 (J=J-N/4),再判 断下一位，依次类推，直到完成高位加1,适2(1 + 1二10b,向右进位 的运算。LH=N/2J=LHI可以从1变换到(N/2-1),Nk=N-2三I J即后半部不用考虑T=X(I) A(IX=X(J) A(J)=TYIN只需前半部和后半部交换后半部不要再重复交换判读各个高位是否为1K = LHJ = J+KJ = J - KK= K/2丄如果该高位为1,则先 减去N/2或N/4、N/8. 再判断下个次高位/输入序列倒序软件程序j = N / 2;/第0个数(二进制数都为0)和最

13、后一个第N4个数(二进制数都为丄)不需倒序 for(i = 1; i N - 2; i+ + )temp = dataRi; dataRi = dataRj; dataRj = temp;k = N / 2;while(l)输入序列倒序的算法设计方法二十进制数X二进制数二进制数十进制数J0000000100110042O1OO1O2W厶2W厶Vni 111 nKJ X XX X KJw67100101110111001 logOil1111537从表格可以看出，所 谓倒序只是把数组下 标的最高位与最低位互换 次高位与次低位互换顺序与倒序二进制数对照表方法二软件分析：已一个字节(N = 256)

14、的倒序为例A0z Alz A255(下标从0000_0000变化到 1111_1111)/*定义两个标志位FO, F1 */ for(i=0;i(255/2);i + + ) /除2是因为只要判断前半部 j=o；F0 = i&0x01;/取最低位Fl=i&0x80;/取最高位if(FO)j=j|0x80; /最低位换到最高位 if(Fl)j=j |0x01; 最高位换到最低位F0=i&0x02;/取次低位Fl=i&0x40;/取次高位if(F0)j=j|0x40; /最次位换到次高位 if(Fl)j=j|0x02; /最次位换到次低位F0 = i&0x04; Fl=i&0x20; if(FO)

15、 j=j|0x20; if(Fl) j=j | 0x04;F0=i&0x08;Fl=i&0xl0; if(FO) j=j|OxlO; if(Fl) j=j|0x08;/前半部与后半部交换，相等时无需交换Ai O AU9只需单层循环即可，共需要循环(128-2)次算法改进一：前面的算法可以进一步优化为:for(i=0;i(255/2);i+ + )for(k=0;k4;k+ + )F0=i&(0x01k);/取最高位if(FO) j=j|(0x80k); /最低位换到最高位 iff(Fl) j=j|(OxOlk); /最高位换到最低位if(ij)/前半部与后半部交换，相等时无需交换Ai O Aj

16、；这个算法只是针对8位的，如果是任意N位的应该如何做? 这里的N必须满足N=2M算法改进二：针对任意N = 2M的情况for(i=0;i(N/2);i + + )/或 #for(i=0;i(N-l)/2);i + + ) j=o；for(k=0;k(M/2);k+ + )/当N=256时，M=8m= 0x01k; FO=i & m; Fl=i & n; if(FO) j=j | n;if(Fl) j=j | m;if(ij)/前半部与后半部交换，相等时无需交换Ai O Aj；F FT软件示例/定义输入波形的幅值#include #define PI 3.1415926 #define N 12

17、8 #define M 7#define AO 255.0void main(void)int izjzk;int p,J,L,B;float SinInN;float dataRN,dataIN; float dataAN;float TrThtemp;/输入波形for(i = 0; i N; i+ + )Sinlni = AO * (sin(2*PI*i/25)+sin(2*PI* i * 0.4 ); dataRi = Sinlni;/输入波形的实数部分datali = 0;/输入波形的虚数部分dataAi = 0;/输出波形的幅度谱为0/输入序列倒序j = N / 2;/第0个数(二进

18、制数都为0)和最后一个第N7个数(二进制数都为丄)不需倒序 for(i = 1; i N - 2; i+ + ) if(i J) temp = dataRi; dataRi = dataRj; dataRU = temp;/因为波形虚数部分都为6所以不用交换/temp = datali;/datali = datalj; /datalO = temp;k = N / 2;while(l) i(j k)/进行FFT/XrJ = Xrf(J) + Tr/XiJ = XQ) + Ti /XrJ+B = Xrf(J) - Tr /XiJ+B = XQ) Ti/(其中X严为上一级的Xiv xr为上一级的

19、Xi)/其中 T= 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) Xrf(J+B)sin(2.0*PI*p/N)for(L=l; L = M;L+ + ) /FFT蝶形级数L从丄一M/计算每个蝶形的两个输入数据相距B = 2(L-1);B = 1;i = L-l; while(i) B = B * 2;/或采用运算，即 B = 2(L-1); B = (int)pow(2,L-l)；/第L级蝶形有pow(2zL-l)r即2的-丄次方个蝶形运算和pow(2zL-l)个旋转因子p for(J=0;J = B-l;J + + )/J =- 1/计算 P = J*2(M-L)p =i = M - L;while(i) p = p * 2;i; p = J * p;/第L级蝶形中同一个旋转因子包含多少个蝶形运算/每个蝶形的两个输入数据相距B = 2(L-1)个点 /同一旋转因子对应着间隔为2八L点的2八(M-L)个蝶形(2八1_= 2*B) for(k = J; k = N-l; k = k+2*B) Tr = da