1、FFT算法簡介php
快速傅里葉變換（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.api
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.
app
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.dom
2、微處理器C語言實現ide
基於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
3、LPC1768硬件仿真lua
1.原始信號3Hz正弦波idea
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三角波spa
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