对于FFT和IFFT的算法和频谱分析的研究

可编辑版/对于FFT和IFFT的算法和频谱分析的研究〔ThealgorithmsandspectrumanalysisofFFTandIFFT摘要：目的在于研究前人的工作结果,对FFT和IFFT有更清楚的认识。主要通过MATLAB的编程完成对FFT和IFFT的算法和频谱分析。首先通过matlab的编程实现FFT和IFFT的这两个函数。然后用已经编译成功的函数实现升余弦滚降。用FFT分析三角函数和三角波函数。用IFFT将上述结果重新变回到时域,通过作图分析变换前后信号的差异。得出了关于fft和ifft函数的分析和关于三角函数和三角波函数的频谱分析的结论关键词：MATLABFFTIFFT升余弦滚降函数三角函数三角波函数Abstract：CompletedthemainalgorithmandspectralanalysisofFFTandIFFTbyMATLABprogramming.First,throughtheMATLABprogrammingtoachievethetwofunctionsFFTandIFFT.Thenusehasbeensuccessfullycompiledfunctionraisedcosine.AnalysisoftrigonometricfunctionandtrianglefunctionbyFFT.WithIFFTtheresultsbackintimedomain,bymappingdifferencesbeforeandaftersignaltransformation.Keywords:MATLAB,FFT,IFFT,Raisedcosinefunction,Trigonometric,Triangularwavefunction引言1965年,库利〔J.W.Cooley和图基〔J.W.Tukey在《计算数学》杂志上发表了"机器计算傅立叶级数的一种算法"的文章,这是一篇关于计算DFT的一种快速有效的计算方法的文章。它的思路建立在对DFT运算内在规律的认识之上。这篇文章的发表使DFT的计算量大大减少,并导致了许多计算方法的发现。这些算法统称为快速傅立叶变换〔FastFourierTransform,简称FFT,1984年,法国的杜哈梅尔〔P.Dohamel和霍尔曼〔H.Hollmann提出的分裂基快速算法,[2]使运算效率进一步提高。FFT即为快速傅氏变换,是离散傅氏变换的快速算法,它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现,但是对于在计算机系统或者说数字系统中应用离散傅立叶变换,可以说是进了一大步。随着科学的进步,FFT算法的重要意义已经远远超过傅里叶分析本身的应用。FFT算法之所以快速,其根本原因在于原始变化矩阵的多余行,此特性也适用于傅里叶变换外的其他一些正交变换,例如,快速沃尔什变换、数论变换等等。在FFT的影响下,人们对于广义的快速正交变换进行了深入研究,使各种快速变换在数字信号处理中占据了重要地位。因此说FFT对数字信号处理技术的发展起了重大推动作用。快速傅里叶变换<FastFourierTranformation,FFT>是将一个大点数N的DFT分解为若干小点的DFT的组合。将运算工作量明显降低,从而大大提高离散傅里叶变换<DFT>的计算速度,从而更加适合进行实时运算。因各个科学技术领域广泛的使用了FFT技术它大大推动了信号处理技术的进步,现已成为数字信号处理强有力的工具,本论文将比较全面的叙述各种快速傅里叶变换算法原理、特点,并完成了基于MATLAB的实现。最后通过FFT和IFFT的两个应用升余弦滚降和确定函数的频谱分析来分别验证FFT和IFFT的正确性和优越性。FFT的算法[1]1.1FFT算法的基本思想设离散的有限长时间序列x<n>,0≤n≤N-1,则其离散傅立叶变换为：这样,矩阵W中有许多相同的元素,从而可以简化DFT的运算过程．FFT算法有许多形式,笔者只讨论最基本的时间抽取基-2FFT算法．1.2算法分析一个N点长序列,直接用DFT方法需要复数乘法N²次；复数加法N<N-1>次。而由图2可知,采用FFT则只需要复数乘法次；复数加法次。当时,这样,运算速度提高了1-2个数量级．图1为FFT算法和直接DFT算法所需运算量与计算点数N的关系曲线．显然,N越大时,优越性越明显．但当N相当大时,利用单机串行进行FFT运算同样满足不了实时系统的需要．[1]1.3算法的程序实现思想及分析 首先检验待变换的序列的元素个数是否为2的幂次方个,如果不是的话则将其补零使之成为2的整数幂次方个。然后利用已经编好的位倒序子程序输出位倒序序号,将输入序列不断分组,进行处理后从新分组,直至完成最后的处理即可输出变换后的结果。需要注意的是,这里fft的变换后结果的元素个数可能与原输入序列的个数不一样,因为如果不是2的幂次方个的话输入序列后面是要补零的。〔主程序为fft_dit和fft_dif1.4流程示意图整个FFT频谱分析与显示过程可用图2所示的流程图示意.fft按频域抽取算法fft按时域抽取算法 图22.IFFT的频谱变换基本原理在实验中为了简化算法我们直接利用前面已经编好的fft_dit或fft_dif这两个现成函数来实现,基本原理如下：将要变换的频谱序列先取共轭然后将其送入前面的函数,将变化后的结果再取共轭即实现IFFT的功能。主程序为ifft_my 图33升3.1升要实现无码间干扰基带传输时,系统必须满足奈奎斯特准则即：对于上述公式,我们分3种情况来说明其含义：Ts<1/2W,其中Ts为系统的输入数据的符号间隔,W为系统的传递函数X〔f的截止频率。由于：因而Z〔f是由频率间隔为1/Ts的X〔f曲线无频率重叠地周期性复制构成。若Ts=1/2W。Z〔f仍是由频率间隔为1/Ts的X〔f曲线无频率重叠地周期性复制构成,在此情况下,仅有一个情况可满足无码间干扰传输的条件,即当此基带传输系统的传递函数是理想低通,其频带宽度为W,则该系统无码间干扰传输的最小Ts=1/2W,即无码间干扰传输的最大符号速率Rs=1/Ts=2W,称此传输速率为奈奎斯特速率。在此理想情况下,虽然系统的频带利用率达到极限,但是此时x<t>是sinc函数,她是非因果的,是物理不可实现的。并且,此x<t>冲击脉冲形状收敛到0的速度极慢,若在收端低通滤波器输出端的采样时科存在定时误差,则在实际采样时刻的采样值会存在码间干扰对于Ts>1/2W情况,Z〔f由频率间隔为1/Ts的X〔f曲线无频率重叠地周期性复制并相加构成的,它还是周期性频谱。在这种情况下,有一特定频谱可满足无码间干扰传输的条件,它就是已获广泛应用的升余弦谱。升余弦滤波器的传递函数表示式为：称α为滚降因子,取值为0≦α≦1。在α=0时,滤波器的带宽W为1/<2Ts>,称为奈奎斯特带宽；α=0.5时,滤波器的截止频率W=〔1+α/〔2Ts=0.75Rs;

α=1时,滤波器的截止频率W=Rs。[3][4]3.2升 升余弦滚降函数是在基带无码间干扰传输中经常用到的频域函数。其主要特性是升余弦滚降函数经过频域平移叠加后能够成一个在各个频域幅度恒定的频域函数。我们在实现升余弦滚降时可以首先在频域实现α=0,α=0.5,α=0.75和α=1四种升余弦滚降函数,然后通过自己编写的ifft_my进行傅里叶逆变换,并将变换后的时域结果反映在图上,分别对比α不同是对应的时域上的不同的时域特性。3.3变换前的时域特性和变换后的频域特性 图44利用fft和ifft进行具体函数频谱分析的实例4.1三角函数的频谱分析及其信号的恢复 在实际信号的分析中,三角函数是非常常见和基本的信号,在这里我们就对三角函数进行分析。在分析中我们会碰到要分析的函数的采样点数不是2的整数次幂个,我们要对它进行补零处理。4.2三角波函数的分析 在分析中因为可能会用到自己编写的T2F函数和F2T函数,但是这两个函数仍然是建立在自己编写的fft_dit和ifft_my的基础上的。只是对输入变量进行了更加完善的处理,一个是补零,另一个是对频域进行了搬移,将原来pi~2pi的部分搬到-pi~0,因为大家在看频谱时比较习惯频谱是关于零对称的。4.3分析结果5结论5.1关于fft和ifft函数的分析 在fft和ifft的实现过程中,的确能降低运算的次数,但是也正好印证了一个很著名的理论,"时间换空间,空间换时间"。在fft和ifft的算法中,我们降低运算的次数是以占用更多的空间换来的。每次将要变换的的序列进行分组,然后对每个组进行处理,虽然降低了运算次数,但是也增加了运算空间的占用。5.2关于升余弦滚降函数的结论 通过图形我们可以观察出,对于不同α的函数,时域主要部分占用的宽度是一样的。但是随着α的增大,时域的起伏是越来越小的。因此我们可以认为,α越大,时域起伏越小,因此在非线性系统中传输时的失真越小。但是与此同时带来的缺点是占用的α越大,频域占用的带宽是越大的,越利于在限带系统中传输,频率的利用率也就越低。5.3关于三角函数和三角波函数的频谱分析结论 〔1根据理论分析,三角函数的频谱因该是几个冲激函数的组合,三角波函数的频谱应该是Sa函数平方的形式。而在实际分析中我们通过观察频谱图可以看出分析的结果基本符合理论分析结果。这也侧面印证了我们自己编写的fft_dit和ifft_my的正确性和有效性。 〔2在通过频谱恢复时域时,我们观察到恢复后的三角函数的时域函数有所失真,三角波函数基本无失真。原因可能是在T2F和F2T函数中可能因为补零改变了元素个数,因此引起了恢复后的函数时域出现失真。致谢：谷群陈若寰王敏附录1：[1]蒋冬初,何飞.FFT算法的并行处理研究.XX城市学院学报<自然科学版>.2005-06-30[2]OppenheimAV,SchaferRW.DigitalSignalProcessing[M].NY:Prentice-Hall,1975.[3]ProakisJG,ManolakisDG.DigitalSignalProcessing-Principles[M].Algorithmsandapplications<2ndEdition>,NY:Printice-Hall,1995.[4]崔灵智.Matlab在数字信号处理课程设计中的应用[J].XX水利职业学院刊,2008.3:11-12.附录2：主要程序和代码升余弦滚降：Fftfunction[num3]=fft_dif<in>N=length<in>;num=in;num3=zeros<1,N>;k=log2<N>;forn=1:kforc=1:2^<n-1>num1=zeros<1,2^<k+1-n>>;num2=zeros<1,2^<k+1-n>>;num1=num<<c-1>*2^<k+1-n>+1:c*2^<k+1-n>>;form=1:2^<k-n>num2<m>=num1<m>+num1<m+2^<k-n>>num2<m+2^<k-n>>=<num1<m>-num1<m+2^<k-n>>>*exp<-i*<m-1>*2*pi/2^<k+1-n>>;endnum<<c-1>*2^<k+1-n>+1:c*2^<k+1-n>>=num2;endendforn=1:Nnum3<n>=num<weidaoxu<n,N>+1>;endfunction[num3]=fft_dit<in>N=length<in>;num=zeros<1,N>;num3=zeros<1,N>;k=log2<N>;forn=1:Nnum<n>=in<weidaoxu<n,N>+1>;endforn=1:kforc=1:2^<k-n>num1=zeros<1,2^n>;num2=zeros<1,2^n>;num1=num<<c-1>*2^n+1:c*2^n>;form=1:2^<n-1>num2<m>=num1<m>+num1<m+2^<n-1>>*exp<-i*<m-1>*2*pi/2^n>;num2<m+2^<n-1>>=num1<m>-num1<m+2^<n-1>>*exp<-i*<m-1>*2*pi/2^n>;endnum<<c-1>*2^n+1:c*2^n>=num2;endend%num3<1>=num<N>;%num3<2:N>=num<1:N-1>;num3=num;IFFTfunction[out]=ifft_my<in>num1=conj<in>;num2=fft_dif<num1>;out=conj<num2>/length<in>;%通过FFT和IFFT来分析升余弦滚降函数的时域和频域的特点clearall;a=0;N=128;l=2;figure<1>[f,out_f]=shengyuxiangunjiang<a,N,l>;subplot<221>;plot<f,out_f>;xlabel<'f'>;ylabel<'a=0时的频谱幅度'>;axis<[-2201.5]>;out_f_=[out_f<<N/2+1>:N>,out_f<1:N/2>];out_t=ifft_my<out_f_>;out_t_=[out_t<<N/2+1>:N>,out_t<1:N/2>];subplot<222>;plot<f,out_t_>;xlabel<'t'>;ylabel<'a=0时的时域幅度'>;axis<[-0.250.25-0.10.5]>;a=0.5;[f,out_f]=shengyuxiangunjiang<a,N,l>;subplot<223>;plot<f,out_f>;xlabel<'f'>;ylabel<'a=0.5时的频谱幅度'>;axis<[-2201.5]>;out_f_=[out_f<<N/2+1>:N>,out_f<1:N/2>];out_t=ifft_my<out_f_>;out_t_=[out_t<<N/2+1>:N>,out_t<1:N/2>];subplot<224>;plot<f,out_t_>;xlabel<'t'>;ylabel<'a=0.5时的时域幅度'>;axis<[-0.250.25-0.10.5]>;a=0.75;figure<2>[f,out_f]=shengyuxiangunjiang<a,N,l>;subplot<221>;plot<f,out_f>;xlabel<'f'>;ylabel<'a=0.75时的频谱幅度'>;axis<[-2201.5]>;out_f_=[out_f<<N/2+1>:N>,out_f<1:N/2>];out_t=ifft_my<out_f_>;out_t_=[out_t<<N/2+1>:N>,out_t<1:N/2>];subplot<222>;plot<f,out_t_>;xlabel<'t'>;ylabel<'a=0.75时的时域幅度'>;axis<[-0.250.25-0.10.5]>;a=1;[f,out_f]=shengyuxiangunjiang<a,N,l>;subplot<223>;plot<f,out_f>;xlabel<'f'>;ylabel<'a=1时的频谱幅度'>;axis<[-2201.5]>;out_f_=[out_f<<N/2+1>:N>,out_f<1:N/2>];out_t=ifft_my<out_f_>;out_t_=[out_t<<N/2+1>:N>,out_t<1:N/2>];subplot<224>;plot<f,out_t_>;xlabel<'t'>;ylabel<'a=1时的时域幅度'>;axis<[-0.250.25-0.10.5]>;逆位倒序：functionout=niweidaoxu<in,long>l=log2<long>;num1=zeros<1,l>;%out=0;fori=1:lnum1<i>=mod<in,2>;out=<in-mod<in,2>>/2;end%fori=1:l%num2<i>=num1<l+1-i>;%end%fori=1:l%out=<2^<i-1>>*num1<i>+out;%endRAND01： functions=rand01<p,m,n> %输入参数: %p:0-1分布中1的概率 %m,n:产生的随机变量样本个数m×n %输出:产生的随机变量样本矢量升余弦滚降： x=rand<m,n>; s=<sign<x-p+eps>+1>/2; %eps=2^<-52>.function[f,out_f]=shengyuxiangunjiang<a,N,l>f=[-l:2*l/N:l-2*l/N];out_f=zeros<1,N>;out_f<1:N/<2*l>*<1-a>>=0;k=[-pi/2:l*pi/<a*N>:pi/2-l*pi/<a*N>]out_f<N/<2*l>*<1-a>+1:N/<2*l>*<1+a>>=<sin<k>+1>/2;out_f<N/<2*l>*<1+a>+1:N/2>=1;form=1:N/2out_f<N+1-m>=out_f<m>;end位倒叙functionout=weidaoxu<in1,long>in=in1-1;l=log2<long>;num1=zeros<1,l>;out=0;fori=1:lnum1<i>=mod<in,2>;in=<in-mod<in,2>>/2;endfori=1:lnum2<i>=num1<l+1-i>;endfori=1:lout=<2^<i-1>>*num2<i>+out;end自余弦滚降function[Ra]=zixiangguan<in>N=length<in>;Ra=zeros<1,N>;form=1:Nfork=1:N-1-mRa<m>=in<k>*in<k+m-1>+Ra<m>;endRa<m>=Ra<m>/<N+1-m>;end频谱分析F变时域function[m]=F2T<M,fs>%输入参数%M：信号的频谱%fs:系统采样频率%输出<返回>参数%m:傅里叶逆变换后的信号,注意其长度为2的整数次幂,利用其画波形时,要注意选取m的一部分,选取长度和所给时间序列t的长度要一致,plot<t,m<1:length<t>>>,否则会出错。m=real<ifft<M>>*fs;FFTfunction[num3]=fft_dif<in1,k>ifnargin==1N=length<in1>;elseN=k;endin=zeros<1,N>;in<1:length<in1>>=in1;num=in;num3=zeros<1,N>;k=log2<N>;forn=1:kforc=1:2^<n-1>num1=zeros<1,2^<k+1-n>>;num2=zeros<1,2^<k+1-n>>;num1=num<<c-1>*2^<k+1-n>+1:c*2^<k+1-n>>;form=1:2^<k-n>num2<m>=num1<m>+num1<m+2^<k-n>>num2<m+2^<k-n>>=<num1<m>-num1<m+2^<k-n>>>*exp<-i*<m-1>*2*pi/2^<k+1-n>>;endnum<<c-1>*2^<k+1-n>+1:c*2^<k+1-n>>=num2;endendforn=1:Nnum3<n>=num<weidaoxu<n,N>+1>;endfunction[num3]=fft_dit<in1,k>ifnargin==1N=length<in1>;elseN=k;endin=zeros<1,N>;in<1:length<in1>>=in1;num=zeros<1,N>;num3=zeros<1,N>;k=log2<N>;forn=1:Nnum<n>=in<weidaoxu<n,N>+1>;endforn=1:kforc=1:2^<k-n>num1=zeros<1,2^n>;num2=zeros<1,2^n>;num1=num<<c-1>*2^n+1:c*2^n>;form=1:2^<n-1>num2<m>=num1<m>+num1<m+2^<n-1>>*exp<-i*<m-1>*2*pi/2^n>;num2<m+2^<n-1>>=num1<m>-num1<m+2^<n-1>>*exp<-i*<m-1>*2*pi/2^n>;endnum<<c-1>*2^n+1:c*2^n>=num2;endend%num3<1>=num<N>;%num3<2:N>=num<1:N-1>;num3=num;function[M,m,df]=fftseq<m,ts,df>%各参数含义与子函数T2F中的完全相同,完成fs=1/ts;ifnargin==2n1=0;elsen1=fs/df;endn2=length<m>;n=2^<max<nextpow2<n1>,nextpow2<n2>>>;M=fft_dit<m,n>;M=[M<n/2+1:n>,M<1:n/2>];m=[m,zeros<1,n-n2>];df=fs/n;function[out]=ifft_my<in>num1=conj<in>;num2=fft_dit<num1>;out=conj<num2>/length<in>;functionout=niweidaoxu<in,long>l=log2<long>;num1=zeros<1,l>;%out=0;fori=1:lnum1<i>=mod<in,2>;out=<in-mod<in,2>>/2;end%fori=1:l%num2<i>=num1<l+1-i>;%end%fori=1:l%out=<2^<i-1>>*num1<i>+out;%endclearalldt=0.001;df=0.2;fs=1/dt;t=[-10:dt:10-dt];y=sin<2*pi*100*t>+2*sin<2*pi*400*t>;[M1,m1,df1,f1]=T2F<y,dt,df,fs>;figure<1>subplot<211>;plot<t,y>;xlabel<'t'>;ylabel<'三角函数频谱分析前的时域'>;axis<[-0.05,0.05,-5,5]>;subplot<212>;plot<f1,abs<M1>>;xlabel<'f'>;ylabel<'三角函数频谱分析后的频谱'>;figure<2>z=zeros<1,length<t>>;z=[zeros<1,length<t>/4>,[0:5*4/length<t>:5-5*4/length<t>