學(xué)術(shù)論文：【畢業(yè)論文】快速傅立葉變換及應(yīng)用

【畢業(yè)論文】快速傅立葉變換及應(yīng)用大連交通大學(xué)2021屆本科生畢業(yè)論文畢業(yè)(設(shè)計)論文題目快速傅里葉變換及其應(yīng)用學(xué)生姓名辛鵬宇專業(yè)班級R計算081班所在院系理學(xué)院指導(dǎo)教師劉立偉職稱副教授所在單位理學(xué)院教研室主任周大勇完成日期2021年6月18日_x000C_摘要快速傅氏變換〔FFT〕，是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn)，但是對于在計算機系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說是進了一大步。傅里葉變換的理論與方法在“數(shù)理方程〞、“線性系統(tǒng)分析〞、“信號處理、仿真〞等很多學(xué)科領(lǐng)域都有著廣泛應(yīng)用,由于計算機只能處理有限長度的離散的序列,所以真正在計算機上運算的是一種離散傅里葉變換。雖然傅里葉運算在各方面計算中有著重要的作用，但是它的計算過于復(fù)雜，大量的計算對于系統(tǒng)的運算負擔過于龐大，使得一些對于耗電量少，運算速度慢的系統(tǒng)對其敬而遠之，然而，快速傅里葉變換的產(chǎn)生，使得傅里葉變換大為簡化，在不犧牲耗電量的條件下提高了系統(tǒng)的運算速度，增強了系統(tǒng)的綜合能力，提高了運算速度，因此快速傅里葉變換在生產(chǎn)和生活中都有著非常重要的作用，對于學(xué)習(xí)掌握都有著非常大的意義。關(guān)鍵字：快速傅氏變換；快速算法；簡化；廣泛應(yīng)用_x000C_ABSTRACTFastFourierTransform(FFT),isadiscretefastFouriertransformalgorithm,whichisbasedontheDiscreteFourierTransformofoddandeven,false,false,andothercharacteristicsoftheDiscreteFourierTransformalgorithmsimprovementsobtained.ItsFouriertransformtheoryhasnotfoundanew,butinthecomputersystemortheapplicationofdigitalsystemsDiscreteFourierTransformcanbesaidtobeabigstepinto.Fouriertransformtheoryandmethodsinthe"mathematicalequation"and"linearsystemsanalysis"and"signalprocessing,simulation,"andmanyotherareashaveawiderangeofapplications,asthecomputercanonlyhandlealimitedlengthofthesequenceofdiscrete,sotrueOnthecomputer'soperationisadiscreteFouriertransform.FourierAlthoughallaspectsofcomputinginthecalculationhasanimportantrole,butitscalculationwastoocomplicated,alotofcomputingsystemforcalculatingtheburdenistoolargeforsomeLesspowerconsumption,theslowspeedofoperationofitssystematarm'slength,however,havethefastFouriertransform,Fouriertransformgreatlysimplifyingthemaking,notinpowerattheexpenseoftheconditionstoincreasethespeedofcomputingsystems,andenhancethesystemThecomprehensiveabilitytoimprovethespeedofoperation,theFastFourierTransformintheproductionandlifehaveaveryimportantroleinlearningtomasterallhavegreatsignificance.KeyWords：FastFourierTransform;fastalgorithm;simplified;widelyused大連交通大學(xué)2021屆本科生畢業(yè)論文18目錄l_Toc31137一、快速傅里葉變換原理及性質(zhì) REF_Toc311371l_Toc28983〔一〕快速傅里葉變換原理 REF_Toc289831l_Toc21714〔二〕快速傅里葉變換的優(yōu)越性 REF_Toc217141l_Toc9754〔三〕快速傅里葉變換的意義 REF_Toc97542l_Toc27753二、快速傅里葉變換的算法 REF_Toc277534l_Toc8124〔一〕快速傅里葉變換算法 REF_Toc81244l_Toc5362三、快速傅里葉變換的應(yīng)用 REF_Toc53626l_Toc25070〔一〕利用FFT計算連續(xù)時間信號的傅里葉變換 REF_Toc250706l_Toc9874〔二〕利用FFT計算離散信號的線性卷積 REF_Toc98749l_Toc401〔三〕利用FFT進行離散信號壓縮 REF_Toc40111l_Toc12997〔四〕利用FFT對離散信號進行濾波 REF_Toc1299714l_Toc25031〔五〕利用FFT提取離散信號中的最強正弦分量 REF_Toc2503117l_Toc12629謝辭 REF_Toc1262922l_Toc32669參考文獻 REF_Toc3266923_x000C_一、快速傅里葉變換原理及性質(zhì)數(shù)字信號的傅里葉變換,通常采用離散傅里葉變換(DFT)方法。DFT存在的缺乏是計算量太大,很難進行實時處理。計算一個N點的DFT,一般需要次復(fù)數(shù)乘法和N(N-1)次復(fù)數(shù)加法運算.因此,當N較大或要求對信號進行實時處理時,往往難以實現(xiàn)所需的運算速度。1965年,J.W.Cooly和J.W.Tukey發(fā)現(xiàn)了DFT的一種快速算法,經(jīng)其他學(xué)者進一步改進,很快形成了一套高效運算方法,這就是現(xiàn)在通用的快速傅里葉變換,簡稱FFT(TheFastFourierTransform)?？焖俑道锶~變換的實質(zhì)是利用式(1)中的權(quán)函數(shù)的對稱性和周期性,把N點DFT進行一系列分解和組合,使整個DFT的計算過程變成一系列疊代運算過程,使DFT的運算量大大簡化,為DFT及數(shù)字信號的實時處理和應(yīng)用創(chuàng)造了良好的條件?！惨弧晨焖俑道锶~變換原理快速傅里葉變換原理1.將長序列DFT分解為短序列的DFT2.利用旋轉(zhuǎn)因子的周期性、對稱性、可約性。將時域序列逐次分解為一組子序列，利用旋轉(zhuǎn)因子的特性，由子序列的DFT來實現(xiàn)整個序列的DFT。其中：快速傅里葉變換分為兩種，分為基2時間抽取算法和基2頻率抽取算法基2時間抽取(Decimationintime)FFT算法其中:r=0,1,2…基2頻率抽取(Decimationinfrequency)FFT算法〔二〕快速傅里葉變換的優(yōu)越性設(shè)為項的復(fù)數(shù)序列，由DFT變換，任一的計算都需要次復(fù)數(shù)乘法和次復(fù)數(shù)加法，而一次復(fù)數(shù)乘法等于四次實數(shù)乘法和兩次實數(shù)加法，一次復(fù)數(shù)加法等于兩次實數(shù)加法，即使把一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法定義成一次“運算〞〔四次實數(shù)乘法和四次實數(shù)加法〕，那么求出項復(fù)數(shù)序列的,即N點DFT變換大約就需要次運算。當點甚至更多的時候，需要次運算，在FFT中，利用WN的周期性和對稱性，把一個N項序列〔設(shè),為正整數(shù)〕，分為兩個項的子序列，每個點DFT變換需要次運算，再用N次運算把兩個點的DFT變換組合成一個N點的DFT變換。這樣變換以后，總的運算次數(shù)就變成。繼續(xù)上面的例子，時，總的運算次數(shù)就變成了525312次，節(jié)省了大約50%的運算量。而如果我們將這種“一分為二〞的思想不斷進行下去，直到分成兩兩一組的DFT運算單元，那么N點的DFT變換就只需要次的運算，在1024點時，運算量僅有10240次，是先前的直接算法的1%，點數(shù)越多，運算量的節(jié)約就越大，這就是FFT的優(yōu)越性.〔三〕快速傅里葉變換的意義傅立葉變換是數(shù)字信號處理領(lǐng)域一種很重要的算法。要知道傅立葉變換算法的意義，首先要了解傅立葉原理的意義。傅立葉原理說明：任何連續(xù)測量的時序或信號，都可以表示為不同頻率的正弦波信號的無限疊加。而根據(jù)該原理創(chuàng)立的傅立葉變換算法利用直接測量到的原始信號，以累加方式來計算該信號中不同正弦波信號的頻率、振幅和相位。和傅立葉變換算法對應(yīng)的是反傅立葉變換算法。該反變換從本質(zhì)上說也是一種累加處理，這樣就可以將單獨改變的正弦波信號轉(zhuǎn)換成一個信號。因此，可以說，傅立葉變換將原來難以處理的時域信號轉(zhuǎn)換成了易于分析的頻域信號〔信號的頻譜〕，可以利用一些工具對這些頻域信號進行處理、加工。最后還可以利用傅立葉反變換將這些頻域信號轉(zhuǎn)換成時域信號。從現(xiàn)代數(shù)學(xué)的眼光來看，傅里葉變換是一種特殊的積分變換。它能將滿足一定條件的某個函數(shù)表示成正弦基函數(shù)的線性組合或者積分。在不同的研究領(lǐng)域，傅里葉變換具有多種不同的變體形式，如連續(xù)傅里葉變換和離散傅里葉變換。在數(shù)學(xué)領(lǐng)域，盡管最初傅立葉分析是作為熱過程的解析分析的工具，但是其思想方法仍然具有典型的復(fù)原論和分析主義的特征。"任意"的函數(shù)通過一定的分解，都能夠表示為正弦函數(shù)的線性組合的形式，而正弦函數(shù)在物理上是被充分研究而相對簡單的函數(shù)類：1.傅立葉變換是線性算子,假設(shè)賦予適當?shù)姆稊?shù),它還是酉算子;2.傅立葉變換的逆變換容易求出,而且形式與正變換非常類似;3.正弦基函數(shù)是微分運算的本征函數(shù),從而使得線性微分方程的求解可以轉(zhuǎn)化為常系數(shù)的代數(shù)方程的求解.在線性時不變雜的卷積運算為簡單的乘積運算,從而提供了計算卷積的一種簡單手段;4.離散形式的傅立葉的物理系統(tǒng)內(nèi),頻率是個不變的性質(zhì),從而系統(tǒng)對于復(fù)雜鼓勵的響應(yīng)可以通過組合其對不同頻率正弦信號的響應(yīng)來獲取;5.著名的卷積定理指出:傅立葉變換可以化復(fù)變換可以利用數(shù)字計算機快速的算出(其算法稱為快速傅立葉變換算法(FFT))。正是由于上述的良好性質(zhì),傅里葉變換在物理學(xué)、數(shù)論、組合數(shù)學(xué)、信號處理、概率、統(tǒng)計、密碼學(xué)、聲學(xué)、光學(xué)等領(lǐng)域都有著廣泛的應(yīng)用。圖像的頻率是表征圖像中灰度變化劇烈程度的指標，是灰度在平面空間上的梯度。如：大面積的沙漠在圖像中是一片灰度變化緩慢的區(qū)域，對應(yīng)的頻率值很低；而對于地表屬性變換劇烈的邊緣區(qū)域在圖像中是一片灰度變化劇烈的區(qū)域，對應(yīng)的頻率值較高。傅立葉變換在實際中有非常明顯的物理意義，設(shè)f是一個能量有限的模擬信號，那么其傅立葉變換就表示f的譜。從純粹的數(shù)學(xué)意義上看，傅立葉變換是將一個函數(shù)轉(zhuǎn)換為一系列周期函數(shù)來處理的。從物理效果看，傅立葉變換是將圖像從空間域轉(zhuǎn)換到頻率域，其逆變換是將圖像從頻率域轉(zhuǎn)換到空間域。換句話說，傅立葉變換的物理意義是將圖像的灰度分布函數(shù)變換為圖像的頻率分布函數(shù)，傅立葉逆變換是將圖像的頻率分布函數(shù)變換為灰度分布函數(shù)。傅立葉變換以前，圖像〔未壓縮的位圖〕是由對在連續(xù)空間〔現(xiàn)實空間〕上的采樣得到一系列點的集合，我們習(xí)慣用一個二維矩陣表示空間上各點，那么圖像可由z=f(x,y)來表示。由于空間是三維的，圖像是二維的，因此空間中物體在另一個維度上的關(guān)系就由梯度來表示，這樣我們可以通過觀察圖像得知物體在三維空間中的對應(yīng)關(guān)系。為什么要提梯度？因為實際上對圖像進行二維傅立葉變換得到頻譜圖，就是圖像梯度的分布圖，當然頻譜圖上的各點與圖像上各點并不存在一一對應(yīng)的關(guān)系，即使在不移頻的情況下也是沒有。傅立葉頻譜圖上我們看到的明暗不一的亮點，實際上圖像上某一點與鄰域點差異的強弱，即梯度的大小，也即該點的頻率的大小〔可以這么理解，圖像中的低頻局部指低梯度的點，高頻局部相反〕。一般來講，梯度大那么該點的亮度強，否那么該點亮度弱。這樣通過觀察傅立葉變換后的頻譜圖，也叫功率圖，我們首先就可以看出，圖像的能量分布，如果頻譜圖中暗的點數(shù)更多，那么實際圖像是比較柔和的〔因為各點與鄰域差異都不大，梯度相對較小〕，反之，如果頻譜圖中亮的點數(shù)多，那么實際圖像一定是鋒利的，邊界清楚且邊界兩邊像素差異較大的。對頻譜移頻到原點以后，可以看出圖像的頻率分布是以原點為圓心，對稱分布的。將頻譜移頻到圓心除了可以清晰地看出圖像頻率分布以外，還有一個好處，它可以別離出有周期性規(guī)律的干擾信號，比方正弦干擾，一副帶有正弦干擾，移頻到原點的頻譜圖上可以看出除了中心以外還存在以某一點為中心，對稱分布的亮點集合，這個集合就是干擾噪音產(chǎn)生的，這時可以很直觀的通過在該位置放置帶阻濾波器消除干擾。二、快速傅里葉變換的算法〔一〕快速傅里葉變換算法快速傅里葉變換算法如下：QUOTE (1)由(1)式可知，對每一個n，計算X(ｎ)須作N次復(fù)數(shù)乘法及N-1次復(fù)數(shù)加法，要完成這組變換共需QUOTE次乘法及N(N-1)次復(fù)數(shù)加法。但以下介紹的快速傅里葉變換的算法，可大大減少運算次數(shù)，提高工作效率。當時，n和k可用二進制數(shù)表示：QUOTE又記,那么〔1〕式可改寫為〔2〕式中：〔3〕因為所以〔2〕可改成〔4〕〔5〕那么式〔5〕即為式〔4〕的分解形式。將初始數(shù)據(jù)代入式(5的第一個等式，可得每一組計算數(shù)據(jù)，一般將痗L-1組計算數(shù)據(jù)代入式(5的第L個等式，計算后可得第L組計算數(shù)據(jù)〔L＝1，2，…，γ〕，計算公式也可表示為=〔6〕式中〔7〕根據(jù)式(5)，第L個數(shù)組中每個的計算只依賴于上一個數(shù)組的兩個數(shù)據(jù)這兩個數(shù)據(jù)的標號相差，即，而且這兩個數(shù)據(jù)只用于計算第L個數(shù)組中標號的數(shù)據(jù)〔等號右端為二進制數(shù)〕。當分別取0和1時，分別有。因此，用上一組的兩個數(shù)據(jù)計算所得的兩個新數(shù)據(jù)仍可儲存在原來位置，計算過程中只需要N個存儲器。將與稱為第L個數(shù)組中的對偶結(jié)點對。計算每個對偶結(jié)點對只需一次乘法，事實上由式〔6〕可得〔8〕式中：；別為式〔7〕中取0，1時對應(yīng)的P值。因，于是對偶結(jié)點的有如下關(guān)系：,因此式〔6〕可表示為〔9〕P的求法：在中，i寫成二進制數(shù)右移位，就成為顛倒位序得式(5)呂，前面的γ個等式，每個等式均對應(yīng)一組數(shù)據(jù)進行計算，每組數(shù)據(jù)都有N/２對結(jié)點，根據(jù)式(9)，每對結(jié)點只需作1次乘法和2次加法，因此，每組數(shù)據(jù)只需N/2次乘法和N次加法，因而完成γ組數(shù)據(jù)的計算共需Nγ/2次乘法和Nγ次加法。三、快速傅里葉變換的應(yīng)用〔一〕利用FFT計算連續(xù)時間信號的傅里葉變換設(shè)是連續(xù)時間信號，并假設(shè)時，那么其傅里葉變換由下式給出令是一個固定的正實數(shù)，是一個固定的正整數(shù)。當時，利用FFT算法可計算。一個固定的時間間隔，選擇足夠小，使得每一個秒的間隔內(nèi)，的變化很小，那么式中積分可近似為〔10〕假設(shè)足夠大，對于所有的整數(shù)，幅值很小，那么式〔10〕變?yōu)椤?1〕當時，式〔10〕兩邊的值為〔12〕其中代表抽樣信號的點。最后令，那么上式變?yōu)椤?3〕首先用FFT算法求出，然后可用上式求出時的。應(yīng)該強調(diào)的是，式〔11〕只是一個近似表示，計算得到的只是一個近似值。通過取更小的抽樣間隔，或者增加點數(shù)，可以得到更精確的值。如果時，幅度譜很小，對應(yīng)于奈奎斯特抽樣頻率，抽樣間隔選擇比較適宜。如果信號只在時間區(qū)間內(nèi)存在，可以通過對時的抽樣信號補零，使足夠大。例1利用FFT計算傅里葉變換如圖1所示的信號其傅里葉變換為：利用下面的命令，可得到的近似值和準確值。圖1連續(xù)時間信號x(t)N=input('InputN:');T=input('InputT:');%計算X(w)近似值t=0:T:2;x=[t-1zeros(1,N-length(t))];X=fft(x);gamma=2*pi/(N*T);k=0:10/gamma;Xapp=(1-exp(-i*k*gamma*T))/(i*k*gamma)*X;%計算真實值X(w)w=0.05:0.05:10;Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w))./(w.*w);plot(k*gamma,abs(Xapp(1:length(k))),'o',w,abs(Xact));legend('近似值','真實值');xlabel('頻率(rad/s)');ylabel('|X|')運行程序后輸入N=128，T=0.1，此時，得到實際的和近似的傅里葉變換的幅度譜如圖2所示，此時近似值已經(jīng)相當準確。通過增加NT可以增加更多的細節(jié)，減少T使得到的值更精確。再次運行程序后輸入N=512，T=0.05，此時，得到實際的和近似的傅里葉變換的幅度譜如圖3所示?！捕忱肍FT計算離散信號的線性卷積兩個離散時間信號與，取，對和右端補零，使得〔14〕利用FFT算法可以求得和的L點DFT，分別是和，利用DTFT卷積性質(zhì)，卷積等于乘積的L點DFT反變換，這也可以通過FFT?算法得到。例2利用FFT計算線性卷積，其中為單位階躍序列，信號如圖4所示。由于當時，很小，故可以取為17；N取10，。利用下面的Matlab命令，可得到、的卷積圖形如圖4所示。subplot(3,1,1);n=0:16;x=0.8.^n;stem(n,x);xlabel('n');ylabel('x[n]');subplot(3,1,2);n=0:15;y=[ones(1,10)zeros(1,6)];stem(n,y);xlabel('n');ylabel('y[n]')subplot(3,1,3);L=26;n=0:L-1;X=fft(x,L);Y=fft(y,L);Z=X.*Y;z=ifft(Z,L);stem(n,z);xlabel('n');ylabel('z[n]')圖4信號x[n]、y[n]及其卷積z[n]=x[n]*y[n]利用下面的Matlab命令，可得到信號x[n]、y[n]的幅度譜與相位譜如圖5所示。subplot(2,2,1);L=26;k=0:L-1;n=0:16;x=0.8.^n;X=fft(x,L);stem(k,abs(X));axis([02505]);xlabel('k');ylabel('|X[k]|')subplot(2,2,2);stem(k,angle(X));axis([025-11]);xlabel('k');ylabel('Angle(X[k])(弧度)')subplot(2,2,3);y=ones(1,10);Y=fft(y,L);stem(k,abs(Y));axis([025010]);xlabel('k');ylabel('|Y[k]|')subplot(2,2,4);stem(k,angle(Y));axis([025-33]);xlabel('k');ylabel('Angle(Y[k])(弧度)')圖5信號x[n]、y[n]的幅度譜與相位譜利用FFT進行離散信號壓縮利用FFT算法對離散信號進行壓縮的步驟如下：1〕通過采樣將信號離散化；2〕對離散化信號進行傅里葉變換；3〕對變換后的系數(shù)進行處理，將絕對值小于某一閾值的系數(shù)置為0，保存剩余的系數(shù)；4〕利用IFFT算法對處理后的信號進行逆傅里葉變換。例3對單位區(qū)間上的以下連續(xù)信號以采樣頻率進行采樣，將其離散化為個采樣值用FFT分解信號，對信號進行小波壓縮，然后重構(gòu)信號。令絕對值最小的80%系數(shù)為0，得到重構(gòu)信號圖形如圖6a)所示，均方差為0.0429，相對誤差為0.0449；令絕對值最小的90%系數(shù)為0，得到重構(gòu)信號圖形如圖6b)所示，均方差為0.0610，相對誤差為0.0638。a)絕對值最小的80%系數(shù)為0的重構(gòu)信號〔FFT〕b)絕對值最小的90%系數(shù)為0的重構(gòu)信號〔FFT〕圖6用FFT壓縮后的重構(gòu)信號相關(guān)Matlab程序如下functionwc=compress(w,r)%輸入信號數(shù)據(jù)w,壓縮率r%輸出壓縮后的信號數(shù)據(jù)if(r<0)|(r>1)error('r應(yīng)該介于0和1之間!');end;N=length(w);Nr=floor(N*r);ww=sort(abs(w));tol=abs(ww(Nr+1));wc=(abs(w)>=tol).*w;function[unbiased_variance,error]=fftcomp(t,y,r)%利用FFT做離散信號壓縮%輸入時間t,原信號y,以及壓縮率r%輸出原信號和壓縮后重構(gòu)信號的圖像,以及重構(gòu)均方差和相對l^2誤差if(r<0)|(r>1)error('r應(yīng)該介于0和1之間!');end;fy=fft(y);yc=ifft(fyc);plot(t,y,'r',t,yc,'b');legend('原信號','重構(gòu)信號');unbiased_variance=norm(y-yc)/sqrt(length(t));error=norm(y-yc)/norm(y);輸入以下Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*sin(8*pi*t);[unbiased_variance,error]=fftcomp(t,f,0.8)unbiased_variance=error=如果用Harr尺度函數(shù)和Harr小波分解信號，對信號進行小波壓縮，然后重構(gòu)信號。令絕對值最小的80%系數(shù)為0，得到重構(gòu)信號圖形如圖7a)所示，均方差為0.0584，相對誤差為0.0611；令絕對值最小的90%系數(shù)為0，得到重構(gòu)信號圖形如圖7b)所示，均方差為0.1136，相對誤差為0.1190。a)絕對值最小的80%系數(shù)為0的重構(gòu)信號〔Harr〕b)絕對值最小的90%系數(shù)為0的重構(gòu)信號〔Harr〕圖7用Harr小波壓縮后的重構(gòu)信號相關(guān)Matlab程序如下function[unbiased_variance,error]=daubcomp(t,y,n,r)%利用Daubechies系列小波做離散信號壓縮%輸入時間t,原信號y,分解層數(shù)n,以及壓縮率r%輸出原信號和壓縮后重構(gòu)信號的圖像,以及重構(gòu)均方差和相對l^2誤差if(r<0)|(r>1)error('r應(yīng)該介于0和1之間!');end;[c,l]=wavedec(y,n,'db1');yc=waverec(cc,l,'db1');plot(t,y,'r',t,yc,'b');legend('原信號','重構(gòu)信號');unbiased_variance=norm(y-yc)/sqrt(length(t));error=norm(y-yc)/norm(y);輸入以下Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*sin(8*pi*t);[unbiased_variance,error]=daubcomp(t,f,8,0.8)unbiased_variance=error=結(jié)論：在信號沒有突變、快變化或者大致上具有周期性的信號，用FFT可以處理得很好〔甚至比小波還要好〕?！菜摹忱肍FT對離散信號進行濾波利用FFT算法對信號進行濾波的步驟如下：1〕通過采樣將信號離散化；2〕對離散化信號進行傅里葉變換；3〕對變換后的系數(shù)進行處理，將絕對值大于某一閾值的系數(shù)置為0，保存剩余的系數(shù)；4〕利用IFFT算法對處理后的信號進行逆傅里葉變換。例4股票價格分析首先進入網(wǎng)址://finance.yahoo/q/hp?s=YHOO，點擊網(wǎng)頁底部位置的DownloadToSpreadsheet按鈕，即可把以Excel表格格式存儲的價格數(shù)據(jù)下載到本地計算機。表格從1列至第6列分別給出了從1996年4月12日至2021年5月30日的交易期里每天的開盤價、最高價、最低價、收盤價、成交量以及趨勢。數(shù)據(jù)下載完成后，需要顛倒順序，使得最早時間的數(shù)據(jù)首先顯示。然后另存到Matlab所在的目錄中，并重新命名為“yhoodata.csv〞。本次分析選擇開盤價，時間是從2021年1月1日至2021年5月30日的=102個交易日期。令代表一支股票的開盤價。為了便于分析，可以先從中減去躍變，得到，即〔15〕輸入以下命令，得到的頻譜如圖8所示。o=csvread('yhoodata.csv',2700,1,[2700128011])';N=102;forn=1:Nx(n)=o(n)-o(1)-((o(N)-o(1))/(N-1))*(n-1);endX=fft(x);k=0:N-1;stem(k,abs(X));axis([01010300]);xlabel('k');ylabel('|X[k]|')圖8x[n]的幅度譜可以看出上圖中有5個較強譜分量，頻率〔〕分別對應(yīng)和。保存這5個頻率分量的系數(shù)，將其他頻率分量的系數(shù)置為0，然后再進行逆傅里葉變換，得到濾波后的近似值。輸入如下Matlab程序，得到真實值與濾波后的近似值，如圖9所示。plot(x);holdon;fliterX=[X(1:2)00X(5)zeros(1,102-9)X(99)00X(102)];fliterx=ifft(fliterX);plot(real(fliterx),'r');axis([110207]);xlabel('n');ylabel('x[n]的值和它的近似值');legend('x[n]真實值','x[n]近似值')圖9x[n]的真實值與濾波后的近似值從上圖可以看出，濾波后的近似值既大致上保存了真實值的變化趨勢，而且與其十分接近。與濾波前比較，濾波后的圖形要比濾波前平滑得多。再由式〔15〕即可求得〔16〕輸入如下Matlab程序，畫出真實開盤價與近似開盤價的圖形。如圖21所示，可以看出是近似根底上的平滑值。plot(o);holdon;forn=1:Noapp(n)=fliterx(n)+o(1)+((o(N)-o(1))/(N-1))*(n-1);endplot(oapp,'r');axis([11022534]);xlabel('n');ylabel('o[n]的值和它的近似值');legend('o[n]真實值','o[n]近似值')圖10o[n]的真實值與濾波后的近似值〔五〕利用FFT提取離散信號中的最強正弦分量這里最強是指在信號中某個正弦分量的振幅遠大于其它正弦分量的振幅?？梢詫η簏c來確定信號中是否有最強的正弦分量。如果信號是連續(xù)時間形式的，首先還要對其進行抽樣，得到離散時間形式的信號，根據(jù)Nyquist定理，抽樣間隔應(yīng)滿，其中是中的最大頻率分量。要判斷信號中是否包含最強正弦分量，采樣數(shù)據(jù)至少要包含該分量一個完整周期的數(shù)據(jù)，例5太陽耀斑數(shù)據(jù)的分析太陽耀斑活動的周期是11年，這個事實可以通過提取耀斑數(shù)據(jù)的最強正弦分量加以證實。耀斑數(shù)據(jù)可以從比利時皇家天文臺(RoyalObservatoryofBelgium)太陽耀斑數(shù)據(jù)索引中心〔SunspotIndexDataCenter,SIDC〕網(wǎng)站下載。網(wǎng)址是：下載后的數(shù)據(jù)存放在文件“monthssn.dat〞中，里面有四列數(shù)據(jù)，第一年是日期，第三列是太陽耀斑的平均數(shù)，第四列平滑后太陽耀斑的平均數(shù)，可以得到從1749年到當前年月〔2021年4月〕的耀斑數(shù)據(jù)。本次分析選擇第三列1981年1月作為開始日期，2021年12月作為結(jié)束日期，共25年300個月份的數(shù)據(jù)。為此先把相關(guān)數(shù)據(jù)復(fù)制到Excel表格的第一列中，然后保存到Matlab所在目錄下，并命名為“sunspotdata.csv〞。然后輸入以下命令，得到耀斑曲線如圖11所示。spd=csvread('sunspotdata.csv',0,0,[002990]);plot(spd);grid;xlabel('月數(shù)');ylabel('耀斑平均數(shù)');axis([03000200])圖111981年1月至2021年12月太陽耀斑的平均數(shù)由上圖可見，太陽耀斑的活動確實具有周期性，但周期的準確值不明顯?？梢酝ㄟ^數(shù)第一個峰值和第二個峰值之間的月份來估計周期的值。查驗表中的數(shù)據(jù)得，第一個峰值為200.3，出現(xiàn)在第116個月〔1990年8月〕，第二個峰值為170.1，出現(xiàn)在第235個月〔2000年7月〕，所以周期是235-116=119月，和實際值132月比較接近。下面利用來分析。首先，從圖11中可以看出從到整個區(qū)間的平均值不為0，為了更方便地分析，需要減去該均值，得到結(jié)果如下：〔17〕然后對進行傅里葉變換，得到。在Matlab中輸入以下命令，得到的幅度譜的圖形如圖12所示。x=spd-sum(spd)/300;X=fft(x);k=0:299;plot(k,abs(X));grid;xlabel('k');ylabel('|X[k]|')圖12的幅度譜由圖12可見，的幅度譜中有一個清晰的尖峰，這就證實中確實包含一個最強正弦分量。為了確定尖峰所對應(yīng)的頻率，用火柴棒圖重畫當時的圖形，在Matlab中輸入以下命令，得到圖形如圖13所示。k=0:10;X=X(1:11);stem(k,abs(X));xlabel('k');ylabel('|X[k]|')圖13當時的火柴棒圖可以看出上圖中有兩個強譜分量，頻率分別為弧度/月和，周期分別為150月和100月。由于的數(shù)據(jù)長度不是其中強正弦分量的整數(shù)倍，譜分量在出現(xiàn)了泄露。要消除泄露，需要使數(shù)據(jù)的長度正好是其中強正弦分量的周期〔太陽耀斑活動的周期，也即〕的整數(shù)倍。為此，重新取中從1至264之間的數(shù)據(jù)進行分析，令〔18〕然后對進行傅里葉變換，得到。輸入以下命令，得到時的火柴棒圖如圖14所示。y=spd(1:264);z=y-sum(y)/264;Z=fft(z);k=0:10;Z=Z(1:11);stem(k,abs(Z));xlabel('k');ylabel('|Z[k]|')圖14當時的火柴棒圖注意到的峰值仍然在，但此時比其附近的值大得多，這說明沒有出現(xiàn)泄露，頻率為弧度/月的正弦分量就是或者中唯一的最強正弦分量。該頻率對應(yīng)周期月，正好等于1