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

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