1、第6講FFT算法設計傅立葉變換將信號從時域轉換為頻域，可以進行 模擬信號的頻率分析離散傅立葉變換（DFT）將信號從頻域轉換為數字 （頻）域，可以進行數字信號（模擬信號數字化） 的頻率分析為了實現DFT在計算機上的快速實現，提出了快速離散傅立葉變換（FFT）如何有傅氏變換->DFT-> FFT?歐拉公式：e*” = cos。土 J sin。e.s.pt.F?) = L /(rk = L feJ dt = )以上5出二-jknX(&) = ZM)e ' ,k =0,1,2,A-lfi=0.2令股:瓦，稱為旋轉_、 N7=> X(k) = 2>()皿如

2、1; = 0,1,2,N-1 =0上式中，k對應數字域，n對應時域另有推導時需用到的公式：T）_/）一泰y J N為I個周期VV N- e- eC-eJ所"'=w； |周期性I二2） w-切-e一片*"'）一一帝"），Nm為加上一個周期-W N-m |周期性|二_ yyN心一聲（1）-聲川-盧也山三3） WN 2 =e N 2 =e J N *e J N 2 =e J N|對稱性=W；"，其中 e ,7T = cos（tt） ysin（7r） = -1.2冗.2萬rX. t 不*4） W， = e n/, =e n =叱1j 可約性推導分

3、析若序列x(n)的長度為N,且滿足N = 2M,(M為自然數) 按n的奇偶性把x(n)分解為兩個N/2的子序列：x1(r)=x(2r)z r=O,T,N/2 - 1x2(r)=x(2r+l)f r=Oflz.vN/2 - 1則x(n)的DFT為：x(k)= Zx()M+=偶數=奇數gN/2T2V/2-1=x(2r)W-kl + £x(2 /+ 1)W產+d r=()r=()/V/2-1N/2-1=£項(廠)£9(/)卬/r=0r=0N/2-1N/2-1Z項w機+w工r=0r=0X(Z) = X(幻+叱江2(，k=0,1,N/2 - 1其中Xx(k)和X2(k)均以

4、N/2為周期NN C NXk + =Xxk + -)+WN 2X2(k-)，k=0zl/-./N/2 - 1NX(k + -)X)-Wx2(k)同理，可推出：Xl(k)=X.(k) + W,.X4(k)n， k=0,l,.”N/4 - 1X + -)=X.(k)-W,2X4(k)X'kj Xsg + WjXJk)X2(k + = X5(k)-W/2X6(k) 9 k=0'，N/4 - 1分到最后，k=0,進行蝶形運算的兩個輸入就是最初輸入序 列x(n)的其中兩個。蝶形分解圖示X(0) X(D X X x(4) x(5) X(6) X(7)N=8點FFT運算圖示x(0) 33x(

5、41A A(3)x(2)雙6)羋機 x(l) X) x(5) 3x(3)上®X _A1Z MD.A(0)A(l)A(2)A(2)A(3)A(4)A(4)AA(5)A(7)A(6)年 X(0) AOI x(l)X 回 X(3) A(4 x(4) 組X 2X(6) 皿XN = T6點FFT運算圖示蝶形運算規律設序列x（n）已經經過時域抽選（倒序）后，存入數組X中。如 果蝶形運算的兩個輸入相距B個點，用原位計算（即計算結果 還放在數組的原來位置），則蝶形運算可表示成如下形式：X,(J) <= X, +.(J + B)W； (x=X,十必X")X/ (J + 8) V=X/

6、|(J) - X/ +(3一 、依十彳)=毛伏)一心x?依)其中：p=J*2ML; J = 0,lf.f2L1-lL=1Z2,.,M下標L表示第L級運算，Xl。）則表示第L級運算后數組元素X。）的值。如果要用實數運算完成上述蝶形運算，可由下面的算法進行:設：T = X,=，+/(3)下標R表示實部X/ (J) = X(/) + jX；(J) (4)下標工表示虛部X%。)代表上一次的實數值Xlt (J + 笈)=X (/ + 6) + jX； (J + B)(5)=> X-(J + 6)W，=Xr(J + 6) + /X/(J + 6)W，= Xr(/ +3) +JX/(/ 4-B)|(c

7、os 手一/sin笥)2tt2tt=XR(J + B)cos v 4 Xj(J + 3)sinp+ jX+ B)cos p X(J + B)sin pNN=Tr + jT, (6)Tr = xr(J + B)cos景 + X； Q + 3)sin素夕 任) q = X/ (J + 8)COS* P-XR(J + 5) sin - PXL(J) = X，J) + jXNJ)= XL_SJ)+TR+jT1由(1)(6)式推導得出=X'R(J) + jX；(J)+TR +jTt由(4)式推導得出XJJ + B) = Xi (/)-(，+ j3)由(2)(6)式推導得出 =Xr (J) + J

8、X/ (/)_( + jTt )由(4)式得出(x/(j)= x；(j)+ r/=> xr(J + B) = Xr(J)-TrXl(J + B) = Xl(J)-T/ (9)公式。)(8)(9)主要用于FFT的軟件編程TH = XH(J + I3)cos27r+ X/(，+ A)sinT, = X,(J + B)cos?- p - XN(J + B)sin2tt72萬一PN送入 x(n), N, M開始輸入序列倒序Xr(J) = Xr(J) + Tr X,(J)=X,(J)+T/Xr(J-B) = Xr(J)-Trx,(J +B) = X/(J)-T/公式中的J就是流程圖中公式的變量k,

9、 流程圖中：N表示階數，M表示總級數，L表示當前級數，B<=2A(L-1)P= J*2A(M-L)X(k) <= X(«)+ X(A+X(k + B) <= Xg - X(k +輸出B表示每個蝶形的兩個輸入數據的間隔, P表示旋轉因子指數可以看出，流程圖總共有3個循環外循環：次數為級數L的變換范圍中循環：為根據當前L求出各個不同的p,循環次 數為P的個數2心工 內循環：為每級中每個P對應的蝶形運算個數, 循環次數為2M內循環中k值每次變換范圍（增量）為2匕這是同一級中每個相同的P對應不同蝶形運算的間隔。看圖推導軟件編程規則：方法1.當N = 8時，第L級共有2工個不

10、同的旋轉因子。因為N = 2M,所以有L=l,2,M,即L的最大值為M（P稱為旋轉因子指數）k=2 （k為p的增量）當 L=1 時,p = 0;當 L=2 時,p = 0,2;當 L=3 時,p = 0,l«23;k=l當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,l,2,3,475,6,7;k=l2 VV （j = 0,&2,3,）=H啜（歸納得出：N/k = 2«J=吧（L= 1,2,3,表當前級數）=WJ"*24一例 （M表總級數）3 .每個p對應每級中的運

11、算個數為:L=4=ML=3L=1L=2有1個蝶形塊，pi=C8個塊 pi=84個蝶形塊有2個蝶形塊pi=4pi=2pi定義為P的增量當L = M時，pi=l = 2° 當乙=例一1時，pi = 2 = 2' 當L = M 2時，pi = 4=22當乙=加=3=1時，pi = 8 = 2反推=>當L = 4H寸，pi = 2當L = 3I時，pi=2M-3 當L = 2H寸，pi = 2"" 當L = 1時，pi = 2"Tp=J*pi=J*2M L （其中J=0,1,2,3,），這樣寫是為了利于軟件的循環編程。此時只要求出每級J的個數J”間

12、即可當L=1時，J= 1 = 2°當L = 2時，兒詢=2 = 2當工=3時,Jtcui = 4 = 22當L = 4時,40tL 8 = 23得：p = J*2ML(J=0,1,2,211)J的總個Tota為2j每一級p的總個數也為2i外循環次數為級數L中循環為根據當前L求出各個不同的p,循環次數為p的個 數27內循環為每級中每個p對應的蝶形運算個數（記為Vtom），循環次數為2Ml當乙=1時，唧0Gd =8 = 23、當L = 2時，VTotaI=4 = 22 =>vTota| = 2m-l當L = 3時，VTotaI=2 = 2'當L = 4時，JTota1 =

13、l = 2° .:每個蝶形的兩個輸入數據間隔（記為INQ ：當L = 1時，W” =1 = 20 ,當L = 2時,1N（I =2 = 2' => |Nd = 21當L = 3時，Wd =4=22|當L = 4時，IN=8 = 23二同一級中每個相同的P對應蝶形運算的間隔（記為VQ ： 當L = 1時，V” = 2 = 2】'當L = 2時，Vz = 4 = 22當L = 3時，匕= 8 = 23=> Vd = 2L當L = 4時，V, =16=24可以看出，為了利于編程，所有的公式推導都往lZ靠輸入序列倒序的算法設計順序倒序十進制數1:二進制數二進制數十

14、進制數J0000000010011004201001023 /I0111106567XVlU 11011 110 111UU JL ion on m1537倒序規律對于N=2M, M位二進制數最高位 的權值為N/2,且從左向右二進制 位的權值依次為N/4,N/8,2,1 o因此，最高位加1相當于十進制運 算J+N/2。如果最高位是0(JvN/2),x(0) x(l) x(2) x(3) x(4) x(5) x(6) x(7)A(0) A(l) 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(

15、4) x(2) x(6) x(l) x(5) x(3) x(7)則直接由J+N/2得下一個倒序值;如果最高位是1(JNN/2),則要將最高位變成0(Jv=JN/2),次高位加1(J+N/4)。但次高位加1時同樣要判斷0、1值，如果為0(JvN/4), 則直接加1(Jv=J+N/4),否則先將次高位變成0(Jv=JN/4),再判斷下一位，依次類推，直到完成高位加1,適2(7 + 7=70。向右進位的運算。LH<=N/2 J<=LH Nk=N-2I可以從1變換到(N/21), 即后半部不用考慮只需前半部和后半部交換 后半部不要再重復交換J <= J - KK<= K/2判讀

16、各個高位是否為1如果該高位為1,則先 減去N/2或N/4、N/8. 再判斷下個次高位/輸入序列倒序軟件程序j = N / 2;/第。個數（二進制數都為0）和最后一個第N-1個數（二進制數都為1）不需倒序for（i = 1; i < N - 2; i+ + ）<temp = dataRi;dataRi = dataRj;dataRj = temp;?k = N / 2;輸入序列倒序的算法設計方法二n而用但1席十進制數工二進制數二進制數十進制數J0000000010011004m nni nni 111 n£KJ X XX X KJ41000011 c367XUl110111

17、Oil 1113 37從表格可以看出，所 謂倒序只是把數組下 標的一一最高位與最低位互換 次高位與次低位互換順序與倒序二進制數對照表方法二軟件分析：己一個字節(N = 256)的倒序為例A0z Alz A255(下標從0000_0000變化至(111工1_工工工1)/*定義兩個標志位FO, F1 */for(i=0;i<(255/2);i + + ) /除2是因為只要判斷前半部 j=0；F0 = i&0x01; 取最低位Fl=i&0x80; 取最高位if(FO)j=j|0x80; /最低位換到最高位if(Fl)j=j|0x01; /最高位換到最低位F0=i&0x0

18、2; 取次低位Fl=i&0x40; 取次高位if(F0)j=j|0x40; /最次位換到次高位if(Fl)j=j|0x02; /最次位換到次低位F0=i&0x04;Fl=i&0x20;if(FO) j=j|0x20;if(Fl) j=j10x04;F0=i&0x08;Fl=i&0xl0;if(FO) j=j|OxlO;if(Fl) j=j|0x08;/前半部與后半部交換，相等時無需交換Ai = Aj;只需單層循環即可，共需要循環(1282)次算法改進一:前面的算法可以進一步優化為:for(i=0;i<(255/2);i+ + )for(k=0;k&

19、lt;4;k+ + )F0=i&(0x01<<k); 取最低位Fl=i&(0x80>>k); 取最高位if(FO) j=j|(0x80»k); /最低位換到最高位iff(Fl) j=j|(OxOl«k); /最高位換到最低位 ?if(i<j)/前半部與后半部交換，相等時無需交換Ai = Aj;這個算法只是針對8位的，如果是任意N位的應該如何做? 這里的N必須滿足N=2M算法改進二：針對任意N = 2M的情況for(i=0;i<(N/2);i + + ) 或者 for(i=0;d(NT)/2);i + + )< J=0

20、；for(k=0;k<(M/2);k+ + )/當N=256時，M=8m= 0x01<<k;n= 0x80>>k;FO=i & m;Fl=i & n;if(FO) j=j | n;if(Fl) j=j | m;)if(i<j)/前半部與后半部交換，相等時無需交換Ai = Aj;F FT軟件示例#include <math.h>#define PI 3.1415926#define N 128#define M 7#define AO 255.0/定義輸入波形的幅值void main(void) <nt Mfk； int P/

21、J/UB; float SinInN; float dataRN,dataIN; float dataAN;float TGTi/temp;/輸入波形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個數(二進制數都為0)和母后一個第NT個數(二進制數都為T)不需倒序for(i = 1; i < N -

22、 2; i+ + ) if(i < j) <temp = dataRi;dataRi = dataRj;dataRU = temp;/因為波形虛數部分都為0,所以不用交換/temp = datali;/datali = datalj; /datalO = temp;)k = N / 2;while(l) < if(J < k)< j = j + k; break;else j = j - k; k = k / 2;?進行FFT/XrJ = Xrf(J) + Tr/XiJ = Xi'。)+ Ti/XrJ+B = Xrf(J) - Tr/XiJ+B = Xi*

23、(J) - Ti/(其中Xr,為上一級的Xr, X為上一級的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)-Xrf(J+B)sin(2.0*PI*p/N)for(L=l; L< = M;L+ + ) /尸4蝶形級數1_從1一乂 </計算每個蝶形的兩個輸入數據相距B = 2(L-1);B = 1;i = L-l;while(i) B = B * 2;；/或采用運算，即 B = 2 人(Ll); B = (int)pow(2,L-l)；第L級蝶形有p。w(2,Ll),即2的L-：L次方個蝶形運算和p。w(2,L-：L)個旋轉因子p for(J=0;J< = B-l;J + + )/J = O,T2,27L-1) - 1</計算 P = J*2 人(M-L)p =；I = M - L;while(i) <p = p * 2; i-; p = J * p;/第L級蝶形中同一個旋轉因子包含多少個蝶形運算