對于FFT和IFFT的算法和頻譜分析的研究

可編輯版/對于FFT和IFFT的算法和頻譜分析的研究〔ThealgorithmsandspectrumanalysisofFFTandIFFT摘要：目的在于研究前人的工作結(jié)果,對FFT和IFFT有更清楚的認識。主要通過MATLAB的編程完成對FFT和IFFT的算法和頻譜分析。首先通過matlab的編程實現(xiàn)FFT和IFFT的這兩個函數(shù)。然后用已經(jīng)編譯成功的函數(shù)實現(xiàn)升余弦滾降。用FFT分析三角函數(shù)和三角波函數(shù)。用IFFT將上述結(jié)果重新變回到時域,通過作圖分析變換前后信號的差異。得出了關于fft和ifft函數(shù)的分析和關于三角函數(shù)和三角波函數(shù)的頻譜分析的結(jié)論關鍵詞：MATLABFFTIFFT升余弦滾降函數(shù)三角函數(shù)三角波函數(shù)Abstract：CompletedthemainalgorithmandspectralanalysisofFFTandIFFTbyMATLABprogramming.First,throughtheMATLABprogrammingtoachievethetwofunctionsFFTandIFFT.Thenusehasbeensuccessfullycompiledfunctionraisedcosine.AnalysisoftrigonometricfunctionandtrianglefunctionbyFFT.WithIFFTtheresultsbackintimedomain,bymappingdifferencesbeforeandaftersignaltransformation.Keywords:MATLAB,FFT,IFFT,Raisedcosinefunction,Trigonometric,Triangularwavefunction引言1965年,庫利〔J.W.Cooley和圖基〔J.W.Tukey在《計算數(shù)學》雜志上發(fā)表了"機器計算傅立葉級數(shù)的一種算法"的文章,這是一篇關于計算DFT的一種快速有效的計算方法的文章。它的思路建立在對DFT運算內(nèi)在規(guī)律的認識之上。這篇文章的發(fā)表使DFT的計算量大大減少,并導致了許多計算方法的發(fā)現(xiàn)。這些算法統(tǒng)稱為快速傅立葉變換〔FastFourierTransform,簡稱FFT,1984年,法國的杜哈梅爾〔P.Dohamel和霍爾曼〔H.Hollmann提出的分裂基快速算法,[2]使運算效率進一步提高。FFT即為快速傅氏變換,是離散傅氏變換的快速算法,它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性,對離散傅立葉變換的算法進行改進獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn),但是對于在計算機系統(tǒng)或者說數(shù)字系統(tǒng)中應用離散傅立葉變換,可以說是進了一大步。隨著科學的進步,FFT算法的重要意義已經(jīng)遠遠超過傅里葉分析本身的應用。FFT算法之所以快速,其根本原因在于原始變化矩陣的多余行,此特性也適用于傅里葉變換外的其他一些正交變換,例如,快速沃爾什變換、數(shù)論變換等等。在FFT的影響下,人們對于廣義的快速正交變換進行了深入研究,使各種快速變換在數(shù)字信號處理中占據(jù)了重要地位。因此說FFT對數(shù)字信號處理技術(shù)的發(fā)展起了重大推動作用?？焖俑道锶~變換<FastFourierTranformation,FFT>是將一個大點數(shù)N的DFT分解為若干小點的DFT的組合。將運算工作量明顯降低,從而大大提高離散傅里葉變換<DFT>的計算速度,從而更加適合進行實時運算。因各個科學技術(shù)領域廣泛的使用了FFT技術(shù)它大大推動了信號處理技術(shù)的進步,現(xiàn)已成為數(shù)字信號處理強有力的工具,本論文將比較全面的敘述各種快速傅里葉變換算法原理、特點,并完成了基于MATLAB的實現(xiàn)。最后通過FFT和IFFT的兩個應用升余弦滾降和確定函數(shù)的頻譜分析來分別驗證FFT和IFFT的正確性和優(yōu)越性。FFT的算法[1]1.1FFT算法的基本思想設離散的有限長時間序列x<n>,0≤n≤N-1,則其離散傅立葉變換為：這樣,矩陣W中有許多相同的元素,從而可以簡化DFT的運算過程．FFT算法有許多形式,筆者只討論最基本的時間抽取基-2FFT算法．1.2算法分析一個N點長序列,直接用DFT方法需要復數(shù)乘法N2次；復數(shù)加法N<N-1>次。而由圖2可知,采用FFT則只需要復數(shù)乘法次；復數(shù)加法次。當時,這樣,運算速度提高了1-2個數(shù)量級．圖1為FFT算法和直接DFT算法所需運算量與計算點數(shù)N的關系曲線．顯然,N越大時,優(yōu)越性越明顯．但當N相當大時,利用單機串行進行FFT運算同樣滿足不了實時系統(tǒng)的需要．[1]1.3算法的程序?qū)崿F(xiàn)思想及分析 首先檢驗待變換的序列的元素個數(shù)是否為2的冪次方個,如果不是的話則將其補零使之成為2的整數(shù)冪次方個。然后利用已經(jīng)編好的位倒序子程序輸出位倒序序號,將輸入序列不斷分組,進行處理后從新分組,直至完成最后的處理即可輸出變換后的結(jié)果。需要注意的是,這里fft的變換后結(jié)果的元素個數(shù)可能與原輸入序列的個數(shù)不一樣,因為如果不是2的冪次方個的話輸入序列后面是要補零的?！仓鞒绦驗閒ft_dit和fft_dif1.4流程示意圖整個FFT頻譜分析與顯示過程可用圖2所示的流程圖示意.fft按頻域抽取算法fft按時域抽取算法 圖22.IFFT的頻譜變換基本原理在實驗中為了簡化算法我們直接利用前面已經(jīng)編好的fft_dit或fft_dif這兩個現(xiàn)成函數(shù)來實現(xiàn),基本原理如下：將要變換的頻譜序列先取共軛然后將其送入前面的函數(shù),將變化后的結(jié)果再取共軛即實現(xiàn)IFFT的功能。主程序為ifft_my 圖33升3.1升要實現(xiàn)無碼間干擾基帶傳輸時,系統(tǒng)必須滿足奈奎斯特準則即：對于上述公式,我們分3種情況來說明其含義：Ts<1/2W,其中Ts為系統(tǒng)的輸入數(shù)據(jù)的符號間隔,W為系統(tǒng)的傳遞函數(shù)X〔f的截止頻率。由于：因而Z〔f是由頻率間隔為1/Ts的X〔f曲線無頻率重疊地周期性復制構(gòu)成。若Ts=1/2W。Z〔f仍是由頻率間隔為1/Ts的X〔f曲線無頻率重疊地周期性復制構(gòu)成,在此情況下,僅有一個情況可滿足無碼間干擾傳輸?shù)臈l件,即當此基帶傳輸系統(tǒng)的傳遞函數(shù)是理想低通,其頻帶寬度為W,則該系統(tǒng)無碼間干擾傳輸?shù)淖钚s=1/2W,即無碼間干擾傳輸?shù)淖畲蠓査俾蔙s=1/Ts=2W,稱此傳輸速率為奈奎斯特速率。在此理想情況下,雖然系統(tǒng)的頻帶利用率達到極限,但是此時x<t>是sinc函數(shù),她是非因果的,是物理不可實現(xiàn)的。并且,此x<t>沖擊脈沖形狀收斂到0的速度極慢,若在收端低通濾波器輸出端的采樣時科存在定時誤差,則在實際采樣時刻的采樣值會存在碼間干擾對于Ts>1/2W情況,Z〔f由頻率間隔為1/Ts的X〔f曲線無頻率重疊地周期性復制并相加構(gòu)成的,它還是周期性頻譜。在這種情況下,有一特定頻譜可滿足無碼間干擾傳輸?shù)臈l件,它就是已獲廣泛應用的升余弦譜。升余弦濾波器的傳遞函數(shù)表示式為：稱α為滾降因子,取值為0≦α≦1。在α=0時,濾波器的帶寬W為1/<2Ts>,稱為奈奎斯特帶寬；α=0.5時,濾波器的截止頻率W=〔1+α/〔2Ts=0.75Rs;

α=1時,濾波器的截止頻率W=Rs。[3][4]3.2升 升余弦滾降函數(shù)是在基帶無碼間干擾傳輸中經(jīng)常用到的頻域函數(shù)。其主要特性是升余弦滾降函數(shù)經(jīng)過頻域平移疊加后能夠成一個在各個頻域幅度恒定的頻域函數(shù)。我們在實現(xiàn)升余弦滾降時可以首先在頻域?qū)崿F(xiàn)α=0,α=0.5,α=0.75和α=1四種升余弦滾降函數(shù),然后通過自己編寫的ifft_my進行傅里葉逆變換,并將變換后的時域結(jié)果反映在圖上,分別對比α不同是對應的時域上的不同的時域特性。3.3變換前的時域特性和變換后的頻域特性 圖44利用fft和ifft進行具體函數(shù)頻譜分析的實例4.1三角函數(shù)的頻譜分析及其信號的恢復 在實際信號的分析中,三角函數(shù)是非常常見和基本的信號,在這里我們就對三角函數(shù)進行分析。在分析中我們會碰到要分析的函數(shù)的采樣點數(shù)不是2的整數(shù)次冪個,我們要對它進行補零處理。4.2三角波函數(shù)的分析 在分析中因為可能會用到自己編寫的T2F函數(shù)和F2T函數(shù),但是這兩個函數(shù)仍然是建立在自己編寫的fft_dit和ifft_my的基礎上的。只是對輸入變量進行了更加完善的處理,一個是補零,另一個是對頻域進行了搬移,將原來pi~2pi的部分搬到-pi~0,因為大家在看頻譜時比較習慣頻譜是關于零對稱的。4.3分析結(jié)果5結(jié)論5.1關于fft和ifft函數(shù)的分析 在fft和ifft的實現(xiàn)過程中,的確能降低運算的次數(shù),但是也正好印證了一個很著名的理論,"時間換空間,空間換時間"。在fft和ifft的算法中,我們降低運算的次數(shù)是以占用更多的空間換來的。每次將要變換的的序列進行分組,然后對每個組進行處理,雖然降低了運算次數(shù),但是也增加了運算空間的占用。5.2關于升余弦滾降函數(shù)的結(jié)論 通過圖形我們可以觀察出,對于不同α的函數(shù),時域主要部分占用的寬度是一樣的。但是隨著α的增大,時域的起伏是越來越小的。因此我們可以認為,α越大,時域起伏越小,因此在非線性系統(tǒng)中傳輸時的失真越小。但是與此同時帶來的缺點是占用的α越大,頻域占用的帶寬是越大的,越利于在限帶系統(tǒng)中傳輸,頻率的利用率也就越低。5.3關于三角函數(shù)和三角波函數(shù)的頻譜分析結(jié)論 〔1根據(jù)理論分析,三角函數(shù)的頻譜因該是幾個沖激函數(shù)的組合,三角波函數(shù)的頻譜應該是Sa函數(shù)平方的形式。而在實際分析中我們通過觀察頻譜圖可以看出分析的結(jié)果基本符合理論分析結(jié)果。這也側(cè)面印證了我們自己編寫的fft_dit和ifft_my的正確性和有效性。 〔2在通過頻譜恢復時域時,我們觀察到恢復后的三角函數(shù)的時域函數(shù)有所失真,三角波函數(shù)基本無失真。原因可能是在T2F和F2T函數(shù)中可能因為補零改變了元素個數(shù),因此引起了恢復后的函數(shù)時域出現(xiàn)失真。致謝：谷群陳若寰王敏附錄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在數(shù)字信號處理課程設計中的應用[J].XX水利職業(yè)學院刊,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來分析升余弦滾降函數(shù)的時域和頻域的特點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> %輸入?yún)?shù): %p:0-1分布中1的概率 %m,n:產(chǎn)生的隨機變量樣本個數(shù)m×n %輸出:產(chǎ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>%輸入?yún)?shù)%M：信號的頻譜%fs:系統(tǒng)采樣頻率%輸出<返回>參數(shù)%m:傅里葉逆變換后的信號,注意其長度為2的整數(shù)次冪,利用其畫波形時,要注意選取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>%各參數(shù)含義與子函數(shù)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<'三角函數(shù)頻譜分析前的時域'>;axis<[-0.05,0.05,-5,5]>;subplot<212>;plot<f1,abs<M1>>;xlabel<'f'>;ylabel<'三角函數(shù)頻譜分析后的頻譜'>;figure<2>z=zeros<1,length<t>>;z=[zeros<1,length<t>/4>,[0:5*4/length<t>:5-5*4/length<t>