一、FFT算法簡介
快速傅里葉變換（Fast Fourier Transform)是離散傅里葉變換的一種快速算法，簡稱FFT，通過FFT可以將一個信號從時域變換到頻域。模擬信號經過A/D轉換變為數字信號的過程稱為采樣。為保證采樣后信號的頻譜形狀不失真，采樣頻率必須大於信號中最高頻率成分的2倍，這稱之為采樣定理。
A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT). Fourier analysis converts a signal from its original domain (often time or space) to a representation in the frequency domain and vice versa. The DFT is obtained by decomposing a sequence of values into components of different frequencies. This operation is useful in many fields, but computing it directly from the definition is often too slow to be practical. An FFT rapidly computes such transformations by factorizing the DFT matrix into a product of sparse (mostly zero) factors. As a result, it manages to reduce the complexity of computing the DFT from {\displaystyle O(n^{2})}
, which arises if one simply applies the definition of DFT, to {\displaystyle O(n\log n)}
, where {\displaystyle n}
is the data size. The difference in speed can be enormous, especially for long data sets where N may be in the thousands or millions. In the presence of round-off error, many FFT algorithms are much more accurate than evaluating the DFT definition directly. There are many different FFT algorithms based on a wide range of published theories, from simple complex-number arithmetic to group theory and number theory.
Fast Fourier transforms are widely used for applications in engineering, science, and mathematics. The basic ideas were popularized in 1965, but some algorithms had been derived as early as 1805. In 1994, Gilbert Strang described the FFT as "the most important numerical algorithm of our lifetime", and it was included in Top 10 Algorithms of 20th Century by the IEEE magazine Computing in Science & Engineering.
The best-known FFT algorithms depend upon the factorization of N, but there are FFTs with O(N log N) complexity for all N, even for prime N. Many FFT algorithms only depend on the fact that {\displaystyle e^{-2\pi i/N}}
is an N-th primitive root of unity, and thus can be applied to analogous transforms over any finite field, such as number-theoretic transforms. Since the inverse DFT is the same as the DFT, but with the opposite sign in the exponent and a 1/N factor, any FFT algorithm can easily be adapted for it.
二、微處理器C語言實現
基於LPC1768最小系統硬件平台，內部模擬產生正弦輸入信號。核心算法C語言部分如下：
#include "..\TKIT_Header\TKIT.h" /* Layer specfication --------------------------------------------------------------------------------------- F(ω) = Fourier[f(t)] = Integral{ f(t)*e^-iwt*dt} 公式描述：公式中F(ω)為f(t)的像函數，f(t)為F(ω)的像原函數。 ------------------------------------------------------------------------------------------------------------- -- Fast Fourier Transform -- and -- Inverse Fast Fourier Transform -- -- 基礎知識點： 信號頻率，F 采用頻率, Fs 采用頻率必須是信號頻率的2倍及以上，才能保證采到的信號沒有失真 -- 采樣獲取到數字信號后，就可以對其做FFT變換了。N個采樣點，經過FFT之后，可以得到N個點的FFT結果，這N個點是以復數 形式存儲的。為了有利於蝶形變換運算，通常N取2的整數次方。 -- FFT后的N點，具有以下幾個物理含義： 1' 第一個點表示0HZ，第N+1個點表示采樣頻率Fs(N+1個點不存在)，從第1個點到N個點，這中間被N-1個點平均分成N等份， 每個點的頻率依次增加。例如某點n所表示的頻率為：Fn=(n-1)*Fs/N 2' FFT能分辨的頻率是：Fs/N，列如，Fs=50，采樣1秒鍾，進行FFT，那么FFT所識別的頻率是1HZ，同樣是Fs=50，采樣兩 秒鍾后進行FFT，那么FFT的分辨率就是0.5HZ. 因此如果要提高頻率分辨力，則必須增加采樣點數，也即采樣時間。頻率分辨率和采樣時間是倒數關系。 3' 由於采樣頻率是數字信號頻率的兩倍及以上，所以我們只需要前N/2個結果即可，從FFT的特性上來看，FFT后N個點是對稱 的，所以也只需要查看前N/2個結果. -- 采樣數據：將ADC采樣的數據按自然序列放在s的實部，虛部為0 -- 假設采樣頻率為fs，采樣點數為N，那么FFT結果就是一個N點的復數，每一個點就對應着一個頻率點， 某一點n(n從1開始)表示的頻率為：fn=(n-1)*fs/N。 舉例說明：用1kHz的采樣頻率采樣128點，則FFT結果的128個數據即對應的頻率點分別是 0,1k/128,2k/128,3k/128,,,,127k/128 Hz -- 這個頻率點的幅值為：該點復數的模值除以N/2（n=1時是直流分量，其幅值是該點的模值除以N）。 假設原始信號的峰值為A，那么FFT的結果的每個點（除了第一個點直流分量之外）的模值就是A的N/2倍。而第一個點就是 直流分量，它的模值就是直流分量的N倍。 -- -- ------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------*/ #if TKIT_FFT_EN //計算復數 #define REALIMG(x) sqrt(FFT_Buffer[x].real*FFT_Buffer[x].real + FFT_Buffer[x].imag*FFT_Buffer[x].imag) int N_FFT=0; //傅里葉變換的點數 int M_of_N_FFT=0; //蝶形運算的級數，N = 2^M int Npart2_of_N_FFT=0; //創建正弦函數表時取PI的1/2 int Npart4_of_N_FFT=0; //創建正弦函數表時取PI的1/4 typedef float ElemType; //原始數據序列的數據類型,可以在這里設置 typedef struct //定義復數結構體 { ElemType real,imag; }TYPE_FFT,*TYPE_FFT_PTR; TYPE_FFT_PTR FFT_Buffer =NULL; ElemType* FFT_BufferSin=NULL; //創建正弦函數表 void CREATE_SIN_TABLE(void) { int i=0; for (i=0; i<=Npart4_of_N_FFT; i++) { FFT_BufferSin[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N); } } ElemType Sin_find(ElemType x) { int i = (int)(N_FFT*x); i = i>>1; if (i>Npart4_of_N_FFT) { i = Npart2_of_N_FFT - i; } return FFT_BufferSin[i]; } ElemType Cos_find(ElemType x) { int i = (int)(N_FFT*x); i = i>>1; if (i<Npart4_of_N_FFT) { return FFT_BufferSin[Npart4_of_N_FFT - i]; } else { return -FFT_BufferSin[i - Npart4_of_N_FFT]; } } //變址 void ExchangeLocation(TYPE_FFT *DataInput) { int nextValue,nextM,i,k,j=0; TYPE_FFT temp; nextValue=N_FFT/2; //變址運算，即把自然順序變成倒位序，采用雷德算法 nextM=N_FFT-1; for (i=0;i<nextM;i++) { if (i<j) //如果i<j,即進行變址 { temp=DataInput[j]; DataInput[j]=DataInput[i]; DataInput[i]=temp; } k=nextValue; //求j的下一個倒位序 while (k<=j) //如果k<=j,表示j的最高位為1 { j=j-k; //把最高位變成0 k=k/2; //k/2，比較次高位，依次類推，逐個比較，直到某個位為0 } j=j+k; //把0改為1 } } //FFT運算函數 void Alg_FFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; ElemType angle; TYPE_FFT W,Temp_XX; ExchangeLocation(FFT_Buffer); for (L=1; L<=M_of_N_FFT; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for (J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; W.imag = -Sin_find(angle); //用C++該函數課聲明為inline W.real = Cos_find(angle); //用C++該函數課聲明為inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for (K=J; K<N_FFT; K=K+step) { KB = K + B; Temp_XX.real = FFT_Buffer[KB].real * W.real-FFT_Buffer[KB].imag*W.imag; Temp_XX.imag = W.imag*FFT_Buffer[KB].real + FFT_Buffer[KB].imag*W.real; FFT_Buffer[KB].real = FFT_Buffer[K].real - Temp_XX.real; FFT_Buffer[KB].imag = FFT_Buffer[K].imag - Temp_XX.imag; FFT_Buffer[K].real = FFT_Buffer[K].real + Temp_XX.real; FFT_Buffer[K].imag = FFT_Buffer[K].imag + Temp_XX.imag; } } } } //IFFT運算函數 void Alg_IFFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; ElemType angle; TYPE_FFT W,Temp_XX; ExchangeLocation(FFT_Buffer); for (L=1; L<=M_of_N_FFT; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for (J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; W.imag = Sin_find(angle); //用C++該函數課聲明為inline W.real = Cos_find(angle); //用C++該函數課聲明為inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for (K=J; K<N_FFT; K=K+step) { KB = K + B; Temp_XX.real = FFT_Buffer[KB].real * W.real-FFT_Buffer[KB].imag*W.imag; Temp_XX.imag = W.imag*FFT_Buffer[KB].real + FFT_Buffer[KB].imag*W.real; FFT_Buffer[KB].real = FFT_Buffer[K].real - Temp_XX.real; FFT_Buffer[KB].imag = FFT_Buffer[K].imag - Temp_XX.imag; FFT_Buffer[K].real = FFT_Buffer[K].real + Temp_XX.real; FFT_Buffer[K].imag = FFT_Buffer[K].imag + Temp_XX.imag; } } } } //初始化FFT程序 //N_FFT是 FFT的點數，必須是2的次方 void Alg_FFTStart(int N_of_FFT) { int i=0; int temp=1; N_FFT = N_of_FFT; //傅里葉變換的點數 ，必須是 2的次方 M_of_N_FFT = 0; //蝶形運算的級數，N = 2^M for (i=0; temp<N_FFT; i++) { temp = 2*temp; M_of_N_FFT++; } Npart2_of_N_FFT = N_FFT/2; //創建正弦函數表時取PI的1/2 Npart4_of_N_FFT = N_FFT/4; //創建正弦函數表時取PI的1/4 FFT_Buffer = (TYPE_FFT_PTR)malloc(N_FFT * sizeof(TYPE_FFT)); FFT_BufferSin = (ElemType *)malloc((Npart4_of_N_FFT+1) * sizeof(ElemType)); CREATE_SIN_TABLE(); //創建正弦函數表 } //結束 FFT運算，釋放相關內存 void Alg_FFTFinish(void) { if (FFT_Buffer != NULL) { free(FFT_Buffer); //釋放內存 FFT_Buffer = NULL; } if (FFT_BufferSin != NULL) { free(FFT_BufferSin); //釋放內存 FFT_BufferSin = NULL; } } //外部輸入采樣信號 void Alg_FFTInput(INT16U i, float real) { if( i>=N_FFT )return; FFT_Buffer[i].real = real; FFT_Buffer[i].imag = 0; } //輸出計算結果 float Alg_FFTGetReal (INT16U i) //讀取實部 { if( i<N_FFT )return FFT_Buffer[i].real; return 0.0; } float Alg_FFTGetImag (INT16U i) //讀取虛部 { if( i<N_FFT )return FFT_Buffer[i].imag; return 0.0; } float Alg_FFTGet (INT16U i) //讀取模值 { if( i<N_FFT )return REALIMG(i); return 0.0; } //產生模擬原始數據輸入 void Alg_FFTInputSimulate(float A,float B,float F, int F_Triangle) { int i,j,k; if( F_Triangle ) { F_Triangle = N_FFT/F_Triangle; //三角波 printf("FFT create input data. real = -%d ~ %d \r\n",F_Triangle/2,F_Triangle/2 ); j = 0; k = 0; for (i=0; i<N_FFT; i++ ,j++,k++)//制造輸入序列 { if( k>F_Triangle )k=0; if( j>=8 ){j=0;Printf("\r\n");} FFT_Buffer[i].real = F_Triangle/2 - k; FFT_Buffer[i].imag = 0; Printf("%5.2uf ",FFT_Buffer[i].real); } Printf("\r\n"); } else { //正弦波 printf("FFT create input data. real = %.2f+%.2f*sin(%.2f*2*PI*i/N_FFT) \r\n",A,B,F); j = 0; for (i=0; i<N_FFT; i++ ,j++)//制造輸入序列 { if( j>=8 ){j=0;Printf("\r\n");} //FFT_Buffer[i].real = sin(PI2*i/N_FFT); FFT_Buffer[i].real = A+B*sin(F*PI2*i/N_FFT); FFT_Buffer[i].imag = 0; Printf("%5.2uf ",FFT_Buffer[i].real); } Printf("\r\n"); } } //主函數 void Alg_FFTDemo(float A,float B,float F, int F_Triangle) { #define FFT_TEST_LEN 64 int i = 0; int j = 0; FP64 *x_array; FP64 *y_array; Printf("FFT demo welcome N=%d!\r\n",FFT_TEST_LEN); x_array = (FP64 *)malloc((FFT_TEST_LEN+1) * sizeof(FP64)); y_array = (FP64 *)malloc((FFT_TEST_LEN+1) * sizeof(FP64)); //初始化各項參數，此函數只需執行一次 Alg_FFTStart(FFT_TEST_LEN); //輸入原始數據 Alg_FFTInputSimulate(A,B,F,F_Triangle); //Input 繪圖 for(i=0;i<N_FFT && i<FFT_TEST_LEN;i++)x_array[i] = FFT_Buffer[i].real; PrintPlot_XY("FFT Input function",TRUE ,x_array,NULL,N_FFT); //進行 FFT計算 Printf("Start FFT:%3dticks.%3dus\r\n",TKIT_TicksGet(),TKIT_TimerCounter()); Alg_FFT(); Printf("Finish FFT:%3dticks.%3dus\r\n",TKIT_TicksGet(),TKIT_TimerCounter()); for (i=0,j=0; i<N_FFT; i++,j++) { if( j>=8 ){j=0;Printf("\r\n");} Printf("%5.2uf ",REALIMG(i)); } Printf("\r\n"); //FFT 繪圖 // for(i=0;i<N_FFT && i<FFT_TEST_LEN;i++)x_array[i] = REALIMG(i); // PrintPlot_XY("Fast Fourier Transform.",TRUE ,x_array,NULL,N_FFT); // for(i=0;i<N_FFT && i<FFT_TEST_LEN;i++) // { // x_array[i] = FFT_Buffer[i].real; // y_array[i] = FFT_Buffer[i].imag; // } // //PrintPlot_XY("Fast Fourier Transform. fx=Real and fy=image",TRUE ,x_array,y_array,N_FFT); //FFT_Buffer[2].real = 0; //FFT_Buffer[2].imag = 0; //進行 IFFT計算 Printf("Start IFFT:%3dticks.%3dus\r\n",TKIT_TicksGet(),TKIT_TimerCounter()); Alg_IFFT(); Printf("Finish IFFT:%3dticks.%3dus\r\n",TKIT_TicksGet(),TKIT_TimerCounter()); for (i=0,j=0; i<N_FFT; i++,j++) { if( j>=8 ){j=0;Printf("\r\n");} Printf("%5.2uf ",FFT_Buffer[i].real/N_FFT); } Printf("\r\n"); //IFFT 繪圖 //for(i=0;i<N_FFT && i<FFT_TEST_LEN;i++)x_array[i] = FFT_Buffer[i].real/N_FFT; //PrintPlot_XY("Inverse Fast Fourier Transform",TRUE ,x_array,NULL,N_FFT); //結束 FFT運算，釋放相關內存 Alg_FFTFinish(); free(x_array); free(y_array); Printf("Close FFT\r\n"); } #endif //TKIT_FFT_EN
三、LPC1768硬件仿真
1.原始信號3Hz正弦波
FFT demo welcome N=64! FFT create input data. real = 0.00+1.00*sin(3.00*2*PI*i/N_FFT) 0.00 0.29 0.55 0.77 0.92 0.99 0.98 0.88 0.70 0.47 0.19 -0.09 -0.38 -0.63 -0.83 -0.95 -1.00 -0.95 -0.83 -0.63 -0.38 -0.09 0.19 0.47 0.70 0.88 0.98 0.99 0.92 0.77 0.55 0.29 0.00 -0.29 -0.55 -0.77 -0.92 -0.99 -0.98 -0.88 -0.70 -0.47 -0.19 0.09 0.38 0.63 0.83 0.95 1.00 0.95 0.83 0.63 0.38 0.09 -0.19 -0.47 -0.70 -0.88 -0.98 -0.99 -0.92 -0.77 -0.55 -0.29
Start FFT:914ticks.323us Finish FFT:918ticks.135us 0.00 0.00 0.00 32.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 32.00 0.00 0.00 Start IFFT:979ticks.340us Finish IFFT:983ticks. 34us -0.00 0.29 0.55 0.77 0.92 0.99 0.98 0.88 0.70 0.47 0.19 -0.09 -0.38 -0.63 -0.83 -0.95 -1.00 -0.95 -0.83 -0.63 -0.38 -0.09 0.19 0.47 0.70 0.88 0.98 0.99 0.92 0.77 0.55 0.29 0.00 -0.29 -0.55 -0.77 -0.92 -0.99 -0.98 -0.88 -0.70 -0.47 -0.19 0.09 0.38 0.63 0.83 0.95 1.00 0.95 0.83 0.63 0.38 0.09 -0.19 -0.47 -0.70 -0.88 -0.98 -0.99 -0.92 -0.77 -0.55 -0.29 Close FFT
1.原始信號3Hz三角波
FFT demo welcome N=64! FFT create input data. real = -10 ~ 10 10.00 9.00 8.00 7.00 6.00 5.00 4.00 3.00 2.00 1.00 0.00 -1.00 -2.00 -3.00 -4.00 -5.00 -6.00 -7.00 -8.00 -9.00 -10.00 -11.00 10.00 9.00 8.00 7.00 6.00 5.00 4.00 3.00 2.00 1.00 0.00 -1.00 -2.00 -3.00 -4.00 -5.00 -6.00 -7.00 -8.00 -9.00 -10.00 -11.00 10.00 9.00 8.00 7.00 6.00 5.00 4.00 3.00 2.00 1.00 0.00 -1.00 -2.00 -3.00 -4.00 -5.00 -6.00 -7.00 -8.00 -9.00
Start FFT:987ticks.415us Finish FFT:991ticks.145us 12.00 21.72 31.67 215.34 20.06 28.69 104.73 19.24 28.86 66.57 17.84 29.42 46.54 16.11 30.04 33.77 14.14 30.61 24.65 11.95 31.10 17.63 9.53 31.48 11.95 6.88 31.77 7.18 3.99 31.94 3.11 1.06 32.00 1.06 3.11 31.94 3.99 7.18 31.77 6.88 11.95 31.48 9.53 17.63 31.10 11.95 24.65 30.61 14.14 33.77 30.04 16.11 46.54 29.42 17.84 66.57 28.86 19.24 104.73 28.69 20.06 215.34 31.67 21.72 Start IFFT:052ticks.511us Finish IFFT:056ticks.103us 10.00 9.00 7.99 6.99 6.00 4.99 3.99 2.99 1.99 0.99 -0.00 -1.00 -1.99 -2.99 -3.99 -5.00 -5.99 -6.99 -7.99 -9.00 -10.00 -11.00 10.00 8.99 8.00 6.99 5.99 5.00 3.99 2.99 1.99 1.00 0.00 -1.00 -1.99 -2.99 -3.99 -5.00 -5.99 -6.99 -8.00 -8.99 -10.00 -11.00 10.00 8.99 7.99 6.99 5.99 4.99 3.99 3.00 1.99 0.99 0.00 -0.99 -1.99 -3.00 -3.99 -5.00 -6.00 -6.99 -7.99 -8.99 Close FFT