C語言的FFT
//----------------------------------------------------------------------------------
//----------------------------------------------------------------------------------
01.void dft()
02.{
03. extern int inv;
04. extern long npt;
05. long k,n;
06. double WN,wk,c,s,XR[size],XI[size];
07. extern complex x[size];
08.
09. WN=2*pi/npt;
10.
11. if(inv==1)
12. WN=-WN;
13.
14. for(k=0;k<npt;++k)
15. {
16. XR[k]=0.0;XI[k]=0.0;
17. wk=k*WN;
18.
19. for(n=0;n<npt;++n)
20. {
21. c=cos(n*wk);s=sin(n*wk);
22. XR[k]=XR[k]+x[n+1].real*c+x[n+1].imag*s;
23. XI[k]=XI[k]-x[n+1].real*s+x[n+1].imag*c;
24. }
25. if(inv==1)
26. {
27. XR[k]=XR[k]/npt;
28. XI[k]=XI[k]/npt;
29. }
30. }
31.
32. for(k=1;k<=npt;++k)
33. {
34. x[k].real=XR[k-1];
35. x[k].imag=XI[k-1];
36. }
37.
38.}
//----------------------------------------------------------------------------------
//----------------------------------------------------------------------------------
離散DFT 的展開式子
圖 (a)為實現這一運算的一般方法,它需要兩次乘法、兩次加減法。考慮到-bW和bW兩個乘法僅相差一負號,可將圖 (a)簡化成圖2.7(b),此時僅需一次乘法、兩次加減法。圖 (b)的運算結構像一蝴蝶通常稱作蝶形運算結構簡稱蝶形結,采用這種表示法,就可以將以上所討論的分解過程用流圖表示。
按照這個辦法,繼續把N/2用2除,由於N=2M,仍然是偶數,可以被2整除,因此可以對兩個N/2點的DFT再分別作進一步的分解。即對{G(k)}和{H(k)}的計算,又可以分別通過計算兩個長度為N/4=2點的DFT,進一步節省計算量,見圖。這樣,一個8點的DFT就可以分解為四個2點的DFT。
如果用反推法,把上圖的蝶形運算FFT,展開求出X(0),其結果和DFT的結果完全一致。