FFT的计算机实现

快速傅立叶变换（FFT）的计算机实现 邓凯 电气07061.FFT运算的原理DFT运算过程中如果有限长序列的长度很长时，即N很大时，在运算过程中所做的乘法和加法运算将很多。即便是采用高速的计算机进行运算，也很难达到信息实时处理的要求。由库利、图基等科学家提出的快速傅里叶变换FFT算法大大减少了运算过程中的乘法和加法次数，适合信息处理对实时性的要求，从而得到广泛应用。按时间抽取的FFT算法的基本思想是将输入的有限长序列首先分成奇数序列和偶数序列，分别计算出奇、偶序列的DFT，然后根据DFT的周期性和对称性质，将其化简，接着将已分成的奇、偶序列再次分别划分成奇、偶序列，即前面的奇序列按其长度再划分成奇数序列和偶数序列；前面的偶序列按其长度再划分成奇数序列和偶数序列。分别计算其DFT后，再按上述方法进行化简。如此反复，直至被划分成的奇、偶序列长度为1为止。1.1基二FFT算法原理若设X(k)=DFTx(n)，且0n, kN-1，N为偶数（如果N为奇数，则添上一个零值点使长度N为偶数）。把它分为奇序列和偶序列： 又 则有 其中， 分别是和的点DFT即要求的值仅仅只需要求 在部分的值即可。这样，我们就可以将一直分为奇序列和偶序列来求其值，直到奇偶序列长度为1为止。1.2蝶形图对于FFT运算，我们通常用蝶形图来表示比如我们表示为下面以N=8为例，运用FFT算法原理计算DFT，如图 1.22、图1.23、图1.24所示。1.3 FFT算法特点（1）原位运算。从图2.9的FFT运算流图中我们可以看出，这种运算是很有规律可循的。其中的每级（每列）运算都是由N/2个蝶形运算所构成，如式（1.31）所示。(式1.31)其中，m表示第m列迭代，k, j表示数据所在的行数，该式所对应的蝶形运算结构如图1.31所示。由图1.31的流程图我们可以看出，对于任一个蝶形中的任何两个节点k和j，经过运算后所得的结果作为下一列（级）k, j两个节点的节点变量，同其他节点变量无关，这样就可以采用原位运算方法。即某一列的N个数据送入存储器经蝶形运算后，结果为另一级（列）数据，而这些数据结果可以以蝶形为单位仍然存储在这同一组存储器中，直到最后输出，中间无需另加存储器。也就是说，蝶形运算的两个输出值仍放回到蝶形的两个输入所在的存储器中（即将原先的输入值覆盖掉）。前一级蝶形运算完成后才进行下一级蝶形原位运算，因而整个运算结构由于采用这种原位运算而大大减少了存储单元的个数，有效地降低了成本。（2）倒位序规律。由图2.9我们还可看出，由于采用了原位计算方法，FFT的输出X(k)是按k值的自然顺序排列在存储单元中的，即排列序列为X(0), X(1), , X(7)，而输入的时间序列都不是按照自然顺序排列，而是按x(0), x(4), x(2), , x(7)，输入到存储单元中。这种排列初看起来像是“混乱无序”的，但实质是有序的。如果我们用二进制数来表示07个数，设为n(n2 n1 n0)2，然后再来观察n (n0 n1 n2)2，如表1.1所示。 十进制数 二进制数二进制数十进制数0000000010011004201001023011110641000011510110156110011371111117由表我们可以看出，只要将 转换成，即将二进制的最高位与最低位相交换、次高位与次低位相交换所得的n就是以自然顺序排列的，故通常n的这种排列规律我们称之为倒位序。 这两条性质给我们编程带来了极大的方便。2C语言实现FFT算法2.1算法基本原理如上所述，FFT算法给编程带来了极大大方便，对比于DFT不仅运算次数大大降低，就连所需的内存空间在运算中也不会增加。主要算法有两部分，一个是倒二进制排序，还有最重要的蝶形算法。倒二进制算法对数组下标进行倒二进制运算，把输入数组按倒二进制排序即可，算法很简单，具体可参见代码蝶形算法核心即为如下公式，只需进行几次循环分级进行蝶形算法即可，具体可参见代码。2.2程序流程图 其中的N表示运算的总级数，即将FFT时域序列分出奇偶序列，直至每个序列点数为1所需的次数。若有1024个点，则N为10。3.FFT的Matlab实验3.1 改变序列的长度N对计算结果的影响x=0.5sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。在Matlab中输入如下代码： clf;fs=100;N=128; %采样频率和数据点数n=0:N-1;t=n/fs; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N); %对信号进行快速Fourier变换mag=abs(y); %求得Fourier变换后的振幅f=n*fs/N; %频率序列subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(N=128);grid on;subplot(2,2,2),plot(f(1:N/2),mag(1:N/2); %绘出Nyquist频率之前随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(N=128);grid on;%对信号采样数据为1024点的处理fs=100;N=1024;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N); %对信号进行快速Fourier变换mag=abs(y); %求取Fourier变换的振幅f=n*fs/N;subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(N=1024);grid on;subplot(2,2,4)plot(f(1:N/2),mag(1:N/2); %绘出Nyquist频率之前随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(N=1024);grid on;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，当 N=128、1024点幅频图如下 fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率01进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致的。为了与真实振幅对应，需要将变换后结果乘以2除以N。由上图可得，点数取得越多，由FFT所得的频谱越精确。4.总结傅里叶变换在物理学、数论、组合数学、信号处理、概率论、统计学、密码学、声学、光学、海洋学、结构动力学等领域都有着广泛的应用（例如在信号处理中，傅里叶变换的典型用途是将信号分解成幅值分量和频率分量）。而FFT则实现了傅里叶变换在计算机领域的广泛应用，学习FFT的原理和实现，对后续课程如电力系统的分析等是很有必要的。本次课程设计，我选择了FFT在计算机上的实现这一个题目，采用C语言编程实现FFT算法，同时用Matlab对C程序进行验证和分析。在完成课设过程中，不仅弄清楚了FFT的原理，还学到了Matlab相关的知识，C语言编程能力也得到了较大的提高。附录 程序清单1.1复数的算法#include complex.hCComplex:CComplex(void)Realpart = 0;ImagPart = 0;CComplex:CComplex(double Real,double Imag)Realpart = Real;ImagPart = Imag;CComplex:CComplex(double real)Realpart = real;ImagPart = 0;CComplex:CComplex(void)/*重载数学运算符 */+重载CComplex CComplex:operator +(CComplex OtherComplex)CComplex temp;temp.Realpart=Realpart+OtherComplex.Realpart;temp.ImagPart=ImagPart+OtherComplex.ImagPart;return temp;/-重载CComplex CComplex:operator -(CComplex OtherComplex)CComplex temp;temp.Realpart=Realpart-OtherComplex.Realpart;temp.ImagPart=ImagPart-OtherComplex.ImagPart;return temp;/*重载CComplex CComplex:operator*(CComplex OtherComplex)CComplex temp;temp.Realpart=Realpart*OtherComplex.Realpart-ImagPart*OtherComplex.ImagPart;temp.ImagPart=Realpart*OtherComplex.ImagPart+ImagPart*OtherComplex.Realpart;return temp;1.2FFT运算#include complex.h#define FFTPOINT 8/*输入时域取样数据*/void InPut(double* pData)double Ipt;for (int i=0 ;iIpt;*(pData+i)=Ipt;/*输出变换后的频域数据，复数形式*/void Output(CComplex* pData)CComplex Out;coutFFT变换后，频域数据为endl;for (int i=0;i=0)printf(%.4fj+%.4fn,Out.ImagPart,Out.Realpart);elseprintf(%.4fj%.4fn,Out.ImagPart,Out.Realpart);/*将输入的数据按倒二进制排序*/void FFTArry(double* pData1,CComplex* pData2)int i=0,m=2,t=0;for (N=0;m = FFTPOINT;N+)m*=2;for(i=0;i=0;j-)if(temp/(int)(pow(2,j)t+=pow(2,N-j-1);temp=temp%(int)(pow(2,j);pData2t.Realpart=*(pData1+i);/extern int N=0;void main()void FFT(CComplex* pData);double FFTDataInFFTPOINT=0.0;CComplex FFTDataOutFFTPOINT;cout请输入原数据,最大个数为：FFTPOINTendl;InPut(FFTDataIn);/输入数据FFTArry(FFTDataIn,FFTDataOut);/将输入数据按到二进制排序FFT(FFTDataOut);/进行FFT变换Output(FFTDataOut);/输出/void FFT(CComplex* pData)int i,r,steplen,k,j;for(i=0;iN;i+) /*实现N级运算*/steplen=(int)(pow(2,i); /*计算步长*/for(r=0;rpow(2,i);r+) CComple