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