第六章傅立葉變換

HarbinInstituteofTechnologyYangkunkyang@1/90第六章快速傅立葉變換6.1離散傅立葉變換及其快速算法6.2FFT計(jì)算流程6.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@2/90第六章快速傅立葉變換

傅立葉變換廣泛應(yīng)用于現(xiàn)代物理中很多領(lǐng)域，包括實(shí)驗(yàn)數(shù)據(jù)的處理，以及理論問(wèn)題的計(jì)算。

一個(gè)物理問(wèn)題可以在“時(shí)間區(qū)域”給與描述，其自變量為t，物理過(guò)程用函數(shù)h(t)表示；也可以在“頻率區(qū)域”上描述，其自變量為f，物理過(guò)程用函數(shù)H(f)表示。

物理學(xué)中把h(t)和H(f)

看作是對(duì)同一物理過(guò)程的兩種不同表示。HarbinInstituteofTechnologyYangkunkyang@3/90第六章快速傅立葉變換

在“時(shí)間區(qū)域”和“頻率區(qū)域”的這兩種表示之間的變換可以通過(guò)傅立葉變換來(lái)表示。HarbinInstituteofTechnologyYangkunkyang@4/90第六章快速傅立葉變換為什么做傅立葉變換？

把一個(gè)時(shí)域的問(wèn)題通過(guò)傅立葉變換轉(zhuǎn)換成頻域的問(wèn)題來(lái)研究，往往會(huì)使問(wèn)題大大簡(jiǎn)化，物理特征更加突出。例如：各種復(fù)雜的振動(dòng)現(xiàn)象，都可以利用傅立葉變換轉(zhuǎn)化為由一些不同頻率、不同振幅的簡(jiǎn)諧振動(dòng)疊加而成。HarbinInstituteofTechnologyYangkunkyang@5/90第六章快速傅立葉變換HarbinInstituteofTechnologyYangkunkyang@6/90第六章快速傅立葉變換從數(shù)學(xué)角度看，傅立葉變換的可分為三類(lèi)方法：解析方法：傅立葉級(jí)數(shù)和傅立葉積分；代數(shù)方法：有限離散傅立葉變換(DFT)；數(shù)值計(jì)算方法：離散傅立葉變換的快速算法(FFT)。DFT:DiscreteFourierTransformFFT:FastFourierTransformHarbinInstituteofTechnologyYangkunkyang@7/906.1離散傅立葉變換及其快速算法6.1.1離散傅立葉變換(DFT)

離散傅立葉變換(DFT)實(shí)際上就是插值方法的一種推廣，就是用eikt作為插值基函數(shù)在插值節(jié)點(diǎn)上進(jìn)行插值。

設(shè)函數(shù)f(t)在0≤t<2π區(qū)間N個(gè)等分節(jié)點(diǎn)2π/N的值為：HarbinInstituteofTechnologyYangkunkyang@8/906.1離散傅立葉變換及其快速算法實(shí)際上，N個(gè)分立點(diǎn)的值常常是采樣的結(jié)果采樣周期：采樣值：HarbinInstituteofTechnologyYangkunkyang@9/906.1離散傅立葉變換及其快速算法插值函數(shù)為：在插值節(jié)點(diǎn)上tn即：用周期為2π的eikt作為插值基函數(shù)HarbinInstituteofTechnologyYangkunkyang@10/906.1離散傅立葉變換及其快速算法用的正交性確定系數(shù)正交性可以表示為HarbinInstituteofTechnologyYangkunkyang@11/476.1離散傅立葉變換及其快速算法當(dāng)當(dāng)可看作首項(xiàng)為1，公比為的等比級(jí)數(shù)正交性證明HarbinInstituteofTechnologyYangkunkyang@12/476.1離散傅立葉變換及其快速算法因此即正交性得以證明。HarbinInstituteofTechnologyYangkunkyang@13/906.1離散傅立葉變換及其快速算法用乘下式，并對(duì)n求和用得HarbinInstituteofTechnologyYangkunkyang@14/906.1離散傅立葉變換及其快速算法上式除皆為零得到：HarbinInstituteofTechnologyYangkunkyang@15/906.1離散傅立葉變換及其快速算法離散傅立葉變換可以表示為：離散傅立葉逆變換為：采樣頻率：頻率間隔：(頻率分辨率)HarbinInstituteofTechnologyYangkunkyang@16/906.1離散傅立葉變換及其快速算法若離散傅立葉變換可以表示為：HarbinInstituteofTechnologyYangkunkyang@17/906.1離散傅立葉變換及其快速算法離散傅立葉變換計(jì)算量1個(gè)操作：一個(gè)復(fù)數(shù)乘法+1個(gè)復(fù)數(shù)加法（求和）1個(gè)節(jié)點(diǎn)F(k)的計(jì)算量：N全部的計(jì)算量：N2HarbinInstituteofTechnologyYangkunkyang@18/906.1離散傅立葉變換及其快速算法離散傅立葉變換計(jì)算量用DFT處理衛(wèi)星照片（1cm2），按1μm分割。處理的節(jié)點(diǎn)數(shù)：計(jì)算量：計(jì)算時(shí)間：每秒億次計(jì)算機(jī)HarbinInstituteofTechnologyYangkunkyang@19/906.1離散傅立葉變換及其快速算法由于DFT計(jì)算量太大，即使依賴(lài)高速計(jì)算機(jī)，也很難實(shí)現(xiàn)對(duì)節(jié)點(diǎn)數(shù)很大的實(shí)際問(wèn)題進(jìn)行實(shí)時(shí)計(jì)算，因此必須對(duì)計(jì)算方法改進(jìn)，減少計(jì)算量。快速傅立葉變換（FFT）：計(jì)算量：計(jì)算時(shí)間：對(duì)上一問(wèn)題計(jì)算量變?yōu)椋篐arbinInstituteofTechnologyYangkunkyang@20/906.1離散傅立葉變換及其快速算法DFT與FFT算法運(yùn)算量比較表NN2N/2×log2N24141648646162563232102480644096192128163844482566553610245122901442304102410485765120HarbinInstituteofTechnologyYangkunkyang@21/906.1離散傅立葉變換及其快速算法6.1.2

快速傅立葉變換計(jì)算過(guò)程(FFT)利用Wkn的下列性質(zhì)節(jié)約計(jì)算量：周期性：對(duì)稱(chēng)性：HarbinInstituteofTechnologyYangkunkyang@22/906.1離散傅立葉變換及其快速算法下面介紹FFT算法設(shè)節(jié)點(diǎn)數(shù)：改寫(xiě)為把DFT這種假設(shè)可以提高算法的效率，且利于編程。HarbinInstituteofTechnologyYangkunkyang@23/906.1離散傅立葉變換及其快速算法為了更加直觀說(shuō)明FFT的思想，取r=3,N=8引入記號(hào)

利用二進(jìn)制數(shù)表示k和n，可以更清楚看出變化規(guī)律。HarbinInstituteofTechnologyYangkunkyang@24/906.1離散傅立葉變換及其快速算法利用周期性對(duì)Wkn簡(jiǎn)化DFT可以用二進(jìn)制寫(xiě)為。HarbinInstituteofTechnologyYangkunkyang@25/906.1離散傅立葉變換及其快速算法周期性HarbinInstituteofTechnologyYangkunkyang@26/906.1離散傅立葉變換及其快速算法傅立葉變換為HarbinInstituteofTechnologyYangkunkyang@27/906.1離散傅立葉變換及其快速算法藍(lán)色部分求和結(jié)果僅依賴(lài)k0，n1，n0，即HarbinInstituteofTechnologyYangkunkyang@28/906.1離散傅立葉變換及其快速算法藍(lán)色部分求和結(jié)果僅依賴(lài)k0，k1，n0，即HarbinInstituteofTechnologyYangkunkyang@29/906.1離散傅立葉變換及其快速算法藍(lán)色部分求和結(jié)果僅依賴(lài)k0，k1，k2，即a3即為傅立葉變換的結(jié)果，但角標(biāo)順序?yàn)槟嫘颉arbinInstituteofTechnologyYangkunkyang@30/906.1離散傅立葉變換及其快速算法可以由下列幾部分完成HarbinInstituteofTechnologyYangkunkyang@31/906.1離散傅立葉變換及其快速算法下面討論計(jì)算量按節(jié)點(diǎn)分解得HarbinInstituteofTechnologyYangkunkyang@32/906.1離散傅立葉變換及其快速算法利用對(duì)稱(chēng)性得HarbinInstituteofTechnologyYangkunkyang@33/906.1離散傅立葉變換及其快速算法從上面計(jì)算過(guò)程可以看出，前四項(xiàng)算出的乘法，可以在后四項(xiàng)的運(yùn)算中使用，即從a0→a1，只需4次復(fù)數(shù)乘法和8次復(fù)數(shù)加法。HarbinInstituteofTechnologyYangkunkyang@34/906.1離散傅立葉變換及其快速算法按節(jié)點(diǎn)分解得HarbinInstituteofTechnologyYangkunkyang@35/906.1離散傅立葉變換及其快速算法利用得從a1→a2，同樣需4次復(fù)數(shù)乘法和8次復(fù)數(shù)加法。HarbinInstituteofTechnologyYangkunkyang@36/906.1離散傅立葉變換及其快速算法按節(jié)點(diǎn)分解得HarbinInstituteofTechnologyYangkunkyang@37/906.1離散傅立葉變換及其快速算法利用得從a2→a3，同樣需4次復(fù)數(shù)乘法和8次復(fù)數(shù)加法。HarbinInstituteofTechnologyYangkunkyang@38/906.1離散傅立葉變換及其快速算法這樣從a0→a3，需復(fù)數(shù)乘法：需復(fù)數(shù)加法：這樣得到N=2r節(jié)點(diǎn)，運(yùn)算量為：HarbinInstituteofTechnologyYangkunkyang@39/906.2FFT計(jì)算流程上一節(jié)介紹的快速傅立葉變換過(guò)程可以用信號(hào)流程圖直觀地表示，對(duì)于N=23節(jié)點(diǎn)的計(jì)算過(guò)程，即從a0→a3計(jì)算過(guò)程，可以用下圖表示：全部計(jì)算層數(shù)為r（N=2r）HarbinInstituteofTechnologyYangkunkyang@40/906.2FFT計(jì)算流程

上圖節(jié)點(diǎn)框內(nèi)的數(shù)值為W的冪指數(shù)，此系數(shù)的確定方法如下：1首先將n用二進(jìn)制表示，如：2將該二進(jìn)制數(shù)右移r-l位，左面空位補(bǔ)零HarbinInstituteofTechnologyYangkunkyang@41/906.2FFT計(jì)算流程3將得到的新二進(jìn)制數(shù)取逆序4再將該二進(jìn)制數(shù)轉(zhuǎn)換為十進(jìn)制，即為W的冪指數(shù)HarbinInstituteofTechnologyYangkunkyang@42/906.2FFT計(jì)算流程HarbinInstituteofTechnologyYangkunkyang@43/906.2FFT計(jì)算流程流程圖特點(diǎn)：所謂對(duì)偶節(jié)點(diǎn)對(duì)，指的是后一列的兩個(gè)節(jié)點(diǎn)僅由前一列的兩個(gè)節(jié)點(diǎn)得到，而與其它節(jié)點(diǎn)無(wú)關(guān)。如：相鄰兩列間，存在對(duì)偶節(jié)點(diǎn)對(duì)。HarbinInstituteofTechnologyYangkunkyang@44/906.2FFT計(jì)算流程對(duì)偶節(jié)點(diǎn)對(duì)的間隔：即如HarbinInstituteofTechnologyYangkunkyang@45/906.2FFT計(jì)算流程若具有N個(gè)節(jié)點(diǎn)，第l層的節(jié)點(diǎn)值為一對(duì)節(jié)點(diǎn)的值，只需進(jìn)行一次乘法和兩次加法運(yùn)算。HarbinInstituteofTechnologyYangkunkyang@46/906.2FFT計(jì)算流程Wp的計(jì)算：令先計(jì)算HarbinInstituteofTechnologyYangkunkyang@47/906.2FFT計(jì)算流程利用三角公式：這樣就得到遞推公式HarbinInstituteofTechnologyYangkunkyang@48/906.2FFT計(jì)算流程練習(xí)構(gòu)照N=16節(jié)點(diǎn)FFT流程圖HarbinInstituteofTechnologyYangkunkyang@49/906.3Matlab離散傅立葉變換函數(shù)fftfft2fftnifftifft2ifftnY=fft(X)Y=fft(X,n)Y=fft(X,n,dim)說(shuō)明:如果X為矩陣，則返回每一列的FFT返回n點(diǎn)FFT，若X的長(zhǎng)度小于n，則末尾補(bǔ)零，大于n，X被刪節(jié)到nDim指定維數(shù)運(yùn)用FFTHarbinInstituteofTechnologyYangkunkyang@50/906.3Matlab離散傅立葉變換函數(shù)連續(xù)信號(hào)的FFT對(duì)于連續(xù)信號(hào)的離散是通過(guò)采樣實(shí)現(xiàn)的。若采樣頻率為理論上則在時(shí)域的采樣時(shí)間間隔為：若時(shí)域的采樣點(diǎn)的數(shù)目為n，則時(shí)域上的節(jié)點(diǎn)坐標(biāo)為：或頻域上的節(jié)點(diǎn)坐標(biāo)為：連續(xù)周期信號(hào)?離散譜連續(xù)非周期信號(hào)?連續(xù)譜HarbinInstituteofTechnologyYangkunkyang@51/906.3Matlab離散傅立葉變換函數(shù)連續(xù)信號(hào)的FFT可能產(chǎn)生的問(wèn)題混疊現(xiàn)象頻譜泄漏柵欄效應(yīng)HarbinInstituteofTechnologyYangkunkyang@52/906.3Matlab離散傅立葉變換函數(shù)混疊現(xiàn)象對(duì)連續(xù)信號(hào)進(jìn)行等間隔采樣時(shí)，如果采樣頻率過(guò)低，采樣后信號(hào)的頻率就會(huì)重疊。這種頻譜的重疊導(dǎo)致的失真稱(chēng)為混疊，而重建出來(lái)的信號(hào)稱(chēng)為原信號(hào)的混疊替身，因?yàn)檫@兩個(gè)信號(hào)有同樣的樣本值。HarbinInstituteofTechnologyYangkunkyang@53/906.3Matlab離散傅立葉變換函數(shù)例6.3.1對(duì)下列信號(hào)采樣分析HarbinInstituteofTechnologyYangkunkyang@54/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@55/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@56/906.3Matlab離散傅立葉變換函數(shù)采樣定理對(duì)于連續(xù)信號(hào)，為保證采樣后能夠不失真地還原信號(hào)，采樣頻率fs必須大于2倍原信號(hào)頻譜的最高截止頻率fc，即例6.3.2對(duì)下列信號(hào)采樣分析HarbinInstituteofTechnologyYangkunkyang@57/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@58/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@59/906.3Matlab離散傅立葉變換函數(shù)對(duì)于頻率為fc的無(wú)限長(zhǎng)的連續(xù)周期信號(hào)，它的頻譜應(yīng)該只是在fc處有離散譜。但是，在利用DFT求它的頻譜時(shí)，對(duì)時(shí)域做了截短，結(jié)果使信號(hào)的頻譜不只是在fc處有離散譜，而是在以fc為中心的頻帶范圍內(nèi)都有譜線出現(xiàn)，它們可以理解為是從fc頻率上“泄露”出去的，這種現(xiàn)象稱(chēng)為頻譜“泄露”。DFT變換頻譜泄漏的根本原因是信號(hào)的截?cái)?。即時(shí)域加窗，對(duì)應(yīng)為頻域卷積，因此，窗函數(shù)的寬度等就會(huì)影響到頻譜。頻譜泄漏HarbinInstituteofTechnologyYangkunkyang@60/906.3Matlab離散傅立葉變換函數(shù)一般在時(shí)域經(jīng)常采用矩形截?cái)?，矩形的傅里葉變換為sinc函數(shù)。這個(gè)過(guò)程在頻域相當(dāng)于頻域函數(shù)與一個(gè)sinc函數(shù)進(jìn)行卷積，卷積的結(jié)果是實(shí)際頻譜的展寬。HarbinInstituteofTechnologyYangkunkyang@61/906.3Matlab離散傅立葉變換函數(shù)頻譜泄漏將導(dǎo)致頻譜分辨率降低，可以通過(guò)增加時(shí)域窗口寬度降低矩形窗的主瓣寬度，提高分辨力。在矩形窗寬度不變的情況下，單純提高采樣頻率無(wú)法提高分辨率。例6.3.3對(duì)下列信號(hào)采樣分析，研究頻譜泄露對(duì)分辨率的影響。HarbinInstituteofTechnologyYangkunkyang@62/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@63/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@64/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@65/906.3Matlab離散傅立葉變換函數(shù)1.無(wú)限長(zhǎng)序列為非周期信號(hào)非周期的無(wú)限長(zhǎng)序列，任意截取一段有限長(zhǎng)的序列，都不能代表實(shí)際信號(hào)，分析結(jié)果當(dāng)然與實(shí)際信號(hào)不一致。2.無(wú)限長(zhǎng)序列為周期信號(hào)周期的無(wú)限長(zhǎng)序列，假設(shè)截取的是正好一個(gè)或整數(shù)個(gè)信號(hào)周期的序列，這個(gè)有限長(zhǎng)序列就可以代表原無(wú)限長(zhǎng)序列，如果分析的方法得當(dāng)?shù)脑?，分析結(jié)果應(yīng)該與實(shí)際信號(hào)一致，可以抑制頻率泄露。從無(wú)限長(zhǎng)序列中截取一個(gè)或整數(shù)個(gè)周期，我們稱(chēng)為整周期截?cái)?，反之，稱(chēng)為非整周期截?cái)?。無(wú)限長(zhǎng)序列信號(hào)可以分為兩種情況：HarbinInstituteofTechnologyYangkunkyang@66/906.3Matlab離散傅立葉變換函數(shù)只有連續(xù)周期信號(hào)截?cái)鄷r(shí)持續(xù)時(shí)間與信號(hào)周期呈整數(shù)倍關(guān)系時(shí)，利用DFT變換可以得到精確的模擬信號(hào)頻譜。HarbinInstituteofTechnologyYangkunkyang@67/906.3Matlab離散傅立葉變換函數(shù)例6.3.4對(duì)下列信號(hào)采樣分析，研究整周期截?cái)嗟挠绊?。HarbinInstituteofTechnologyYangkunkyang@68/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@69/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@70/906.3Matlab離散傅立葉變換函數(shù)減少頻率泄漏的方法（1）適當(dāng)加大窗口寬度（2）采用適當(dāng)形狀的加權(quán)窗函數(shù)，加權(quán)的窗函數(shù)包括平頂窗、漢寧窗、高斯窗等等。理想的窗函數(shù)是傅里葉變換后主瓣很窄，旁瓣衰減很快。矩形窗的主瓣很窄，但是旁瓣衰減卻很慢，hanning窗、hamming窗、blackman窗等的旁瓣衰減有了明顯的改進(jìn)，但是主瓣卻寬了很多，大概是矩形窗主瓣的二倍，blackman窗的主瓣還要寬，這就造成了信號(hào)頻譜的頻率識(shí)別率很低！HarbinInstituteofTechnologyYangkunkyang@71/906.3Matlab離散傅立葉變換函數(shù)利用DFT逼近連續(xù)時(shí)間信號(hào)的傅立葉變換,其頻譜將不再是連續(xù)函數(shù)而是基頻f的整數(shù)倍。用DFT計(jì)算頻譜,就如通過(guò)一個(gè)柵欄觀看一個(gè)景色,只能在離散點(diǎn)的地方看到真實(shí)的景象,從而產(chǎn)生柵欄效應(yīng)。

如果在兩離散的譜線間頻譜有很大變化,不作特殊處理,則無(wú)法將其檢測(cè)出來(lái)。柵欄效應(yīng)HarbinInstituteofTechnologyYangkunkyang@72/906.3Matlab離散傅立葉變換函數(shù)減小柵欄效應(yīng)的一個(gè)方法是在所取數(shù)據(jù)的末端加一些零值點(diǎn),使一個(gè)周期內(nèi)點(diǎn)數(shù)增加,但是不改變?cè)械挠涗洈?shù)據(jù)。

這種方法等效于減小了頻率間隔，因頻率間隔為fs/N(fs為采樣頻率)。

增加零值點(diǎn)后,N增加,從而能保持原來(lái)頻譜形式不變的情況下使譜線變密,也就使頻譜抽樣點(diǎn)數(shù)增加。

這樣,原來(lái)看不到的頻譜分量就有可能看到了。補(bǔ)加零點(diǎn)以改變周期時(shí),所用窗函數(shù)寬度卻不能變,亦即必須按數(shù)據(jù)記錄原長(zhǎng)來(lái)選擇窗函數(shù),而不能按補(bǔ)了零值點(diǎn)后的長(zhǎng)度來(lái)選擇窗函數(shù)。

通俗地說(shuō),就是應(yīng)先加窗,再補(bǔ)零。減小柵欄效應(yīng)方法HarbinInstituteofTechnologyYangkunkyang@73/906.3Matlab離散傅立葉變換函數(shù)例6.3.5對(duì)下列信號(hào)采樣分析，研究補(bǔ)加零點(diǎn)的影響。HarbinInstituteofTechnologyYangkunkyang@74/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@75/906.3Matlab離散傅立葉變換函數(shù)HarbinInstituteofTechnologyYangkunkyang@76/906.3Matlab離散傅立葉變換函數(shù)傅立葉變換的正頻率和負(fù)頻率傅立葉變換是以復(fù)指數(shù)eiωt為基函數(shù)展開(kāi)的，實(shí)函數(shù)可以展開(kāi)為兩個(gè)復(fù)指數(shù)信號(hào)的和。的傅立葉變換HarbinInstituteofTechnologyYangkunkyang@77/906.3Matlab離散傅立葉變換函數(shù)所以，對(duì)于實(shí)函數(shù)傅立葉變換后產(chǎn)生兩個(gè)頻率正頻率：ω

負(fù)頻率：-ω

HarbinInstituteofTechnologyYangkunkyang@78/906.3Matlab離散傅立葉變換函數(shù)1）實(shí)信號(hào)FFT的結(jié)果前半部分[0,fs/2]對(duì)應(yīng)的是正頻率，后半部分[-fs/2,0]對(duì)應(yīng)的是負(fù)頻率的結(jié)果。大于fs/2的部分的頻譜實(shí)際上是實(shí)信號(hào)的負(fù)頻率加fs的結(jié)果。2）如果要讓實(shí)信號(hào)FFT的結(jié)果與[-fs/2,fs/2]對(duì)應(yīng)，則需要利用函數(shù)fftshift處理一下即可。fftshift(X):將零頻點(diǎn)移到頻譜的中間ifftshift(X)HarbinInstituteofTechnologyYangkunkyang@79/906.3Matlab離散傅立葉變換函數(shù)例6.3.6實(shí)信號(hào)FFT變換處理。HarbinInstituteofTechnologyYangkunkyang@80/906.3Matla