第六章傅立葉變換

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

傅立葉變換廣泛應用于現(xiàn)代物理中很多領域，包括實驗數據的處理，以及理論問題的計算。

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

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

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

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

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

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

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

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

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

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

如果在兩離散的譜線間頻譜有很大變化,不作特殊處理,則無法將其檢測出來。柵欄效應HarbinInstituteofTechnologyYangkunkyang@72/906.3Matlab離散傅立葉變換函數減小柵欄效應的一個方法是在所取數據的末端加一些零值點,使一個周期內點數增加,但是不改變原有的記錄數據。

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

增加零值點后,N增加,從而能保持原來頻譜形式不變的情況下使譜線變密,也就使頻譜抽樣點數增加。

這樣,原來看不到的頻譜分量就有可能看到了。補加零點以改變周期時,所用窗函數寬度卻不能變,亦即必須按數據記錄原長來選擇窗函數,而不能按補了零值點后的長度來選擇窗函數。

通俗地說,就是應先加窗,再補零。減小柵欄效應方法HarbinInstituteofTechnologyYangkunkyang@73/906.3Matlab離散傅立葉變換函數例6.3.5對下列信號采樣分析，研究補加零點的影響。HarbinInstituteofTechnologyYangkunkyang@74/906.3Matlab離散傅立葉變換函數HarbinInstituteofTechnologyYangkunkyang@75/906.3Matlab離散傅立葉變換函數HarbinInstituteofTechnologyYangkunkyang@76/906.3Matlab離散傅立葉變換函數傅立葉變換的正頻率和負頻率傅立葉變換是以復指數eiωt為基函數展開的，實函數可以展開為兩個復指數信號的和。的傅立葉變換HarbinInstituteofTechnologyYangkunkyang@77/906.3Matlab離散傅立葉變換函數所以，對于實函數傅立葉變換后產生兩個頻率正頻率：ω

負頻率：-ω

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