04快速傅里葉變換

1、為了分析連續(xù)信號頻譜，需要使用為了分析連續(xù)信號頻譜，需要使用dft來逼近連續(xù)時間信來逼近連續(xù)時間信號的傅里葉變換。但由于號的傅里葉變換。但由于dft要對連續(xù)信號進行采樣和截要對連續(xù)信號進行采樣和截斷等處理，因此，常產(chǎn)生誤差分析，體現(xiàn)在三個方面。斷等處理，因此，常產(chǎn)生誤差分析，體現(xiàn)在三個方面。 混疊現(xiàn)象：混疊現(xiàn)象： 截斷效應(yīng)：截斷效應(yīng)： 柵欄效應(yīng)：柵欄效應(yīng)：時域采樣不滿足奈奎斯特準則：時域采樣不滿足奈奎斯特準則： 其中其中 為抽樣頻率，為抽樣頻率， 為信號最高頻率。為信號最高頻率。此條件只規(guī)定了此條件只規(guī)定了 的下限，其上限要受的下限，其上限要受抽樣抽樣間隔）間隔） 的約束。的約束。 也可表示為

2、記錄長度的倒數(shù)，即：也可表示為記錄長度的倒數(shù)，即：其中，其中， 為抽樣點數(shù)。為抽樣點數(shù)。021shpffftnnsfhfsfn2shff0f0f采樣前對信號進行預(yù)濾波，采樣時滿足采樣定理，以采樣前對信號進行預(yù)濾波，采樣時滿足采樣定理，以避免信號在避免信號在w= 處附近的混迭，再濾去信號中頻率高于處附近的混迭，再濾去信號中頻率高于的頻率分量。的頻率分量。2sf 當(dāng)當(dāng) 給定時給定時, 若為了避免混疊現(xiàn)象而一味提高抽樣頻率若為了避免混疊現(xiàn)象而一味提高抽樣頻率 , 導(dǎo)致導(dǎo)致 增加增加，即導(dǎo)致即導(dǎo)致頻率分辨能力下降頻率分辨能力下降; 反之反之, 若打算要提高頻率分辨率即減小若打算要提高頻率分辨率即減小

3、, 則導(dǎo)致減小則導(dǎo)致減小 , 最終必須減最終必須減小信號的高頻容量。小信號的高頻容量。顯然，針對高頻容量顯然，針對高頻容量 與頻率分辨率與頻率分辨率 , 想保持其中一個不變而使想保持其中一個不變而使另一個性能得以提高的唯一辦法是：另一個性能得以提高的唯一辦法是：增加記錄長度內(nèi)的點數(shù)增加記錄長度內(nèi)的點數(shù)n, 即：即：nsf0f0fsfhf0f02hfnf( )x n( )( ) ( )nxnx n w n( )w n()jwx e()jww e()()*()jwjwjwnxex ew e窗函數(shù)不可能取無限寬，即其頻譜不可能為一沖激函數(shù)，窗函數(shù)不可能取無限寬，即其頻譜不可能為一沖激函數(shù)，信號的頻譜

4、與窗函數(shù)的卷積信號的頻譜與窗函數(shù)的卷積如下圖所示。如下圖所示。002/sa， 求極限：的連續(xù)頻譜。為處的一根譜線，變成了以原來在為中心，形狀 1lim/2lim221lim/22jjjjsax ewjx esax ex ejwjw 的窗函數(shù)。接近為了盡量減少泄漏，需要尋找頻域中窗函數(shù)具有旁瓣，并且主瓣也不過，泄漏的產(chǎn)生是由于所以不能用無限加寬窗口來減少泄漏。，意味著需要無限加寬窗函數(shù)，等于未截短。，意味著需要無限加寬窗函數(shù)，等于未截短。若使泄漏若使泄漏至至0，則，則占有一定寬度。，即)(jewnnn2n4n2主瓣旁瓣旁瓣0n點點dft是在區(qū)間是在區(qū)間0, 2 上的上的n點等間隔采樣，點等間隔采

5、樣，就好像從，就好像從n+1個柵欄縫隙中觀看信號的頻譜特性，得到個柵欄縫隙中觀看信號的頻譜特性，得到的是的是n個縫隙中看到的頻譜函數(shù)值，這種現(xiàn)象稱為柵欄效應(yīng)。個縫隙中看到的頻譜函數(shù)值，這種現(xiàn)象稱為柵欄效應(yīng)。對信號的頻譜進行有限點采樣；對信號的頻譜進行有限點采樣；柵欄效應(yīng)可能漏掉（擋住）大的頻譜分量；柵欄效應(yīng)可能漏掉（擋住）大的頻譜分量；對原序列補對原序列補0，增大，增大n，以增加采樣點。，以增加采樣點。減小柵欄效應(yīng)的一個方法是，減小柵欄效應(yīng)的一個方法是，使，使, 但是但是。這種方法這種方法加長了周期加長了周期t0，t0增加增加, 使得使得, 從而能保持原來從而能保持原來的情況下的情況下, 也就

6、使頻譜抽樣點數(shù)增加。也就使頻譜抽樣點數(shù)增加。這樣，原來看不到的頻譜分量就有可能看到了。這樣，原來看不到的頻譜分量就有可能看到了。2f128ssffff88snff( )x n81116542 30109712f2后面，則有：1212snff 2213-j-j400222173-j-j-j88111000( )( )( )10, 1j, 0, 1j( )( )( )( )10, 2.7-j6.5, 1j, 1.3-j0.5, 0, 1.3j0.5, 1j, 2.7j6.5nnmnmnnnnnknknknnnnx kx n ex n ex kx n ex n ex n e 1( )2, 3, 3,

7、 2( )2, 3, 3, 2, 0, 0, 0, 0 x nx n補加零點以改變周期時補加零點以改變周期時, 。, 而不能按補了零值點而不能按補了零值點后的長度來選擇窗函數(shù)。后的長度來選擇窗函數(shù)。即：即：“先加窗先加窗, , 再補零再補零?！币话銇碚f，信號長度 越長，抽樣點數(shù) 越大，則 越小，即分辨力越強。但此處的 是指實際的信號長度， 是指這個長度上的抽樣點數(shù)，而。011=spffnntt采樣頻率采樣點數(shù) 所謂頻率分辨率，是指能夠?qū)⑿盘栔袃蓚€靠得很近的譜峰分開的能力：ptn0fptnt其中：時域采樣間隔1111psttnf時域采樣間隔：11 spnft，所以2221psttnf22spnf

8、t11011122022211sppsppfnfnt ntfnfnt ntpt1sf2sf當(dāng)采樣頻率由變?yōu)闀r，在數(shù)據(jù)長度1n2n變?yōu)?，。不變情況下，采樣點數(shù)由僅增加采樣率并不能改變頻率分辨率。如果數(shù)據(jù)長度不變，則有：pt0102ff已知某信號如下：假設(shè)最高頻率 ，請用 對其抽樣，并觀察其頻譜情況。 假設(shè)取信號長度為：假設(shè)取信號長度為： ，經(jīng)抽樣所得的，經(jīng)抽樣所得的 的點數(shù)為：的點數(shù)為：對對 做做dft變換，所得到的頻率最大分辨率為：變換，所得到的頻率最大分辨率為：( )sin(4)sin(4.04)sin(4.14)x tttt 3hf 10sf 25.6 pts 256psntf ( )x

9、n( )x n01010.0390625 hz256spffnt 只顯示了兩個頻率：f1、f3 f2-f1=0.02f0，所以能分辨出來 增加采樣點數(shù) n，即增加數(shù)據(jù)長度。設(shè)原數(shù)據(jù)長度設(shè)原數(shù)據(jù)長度t1，抽樣點數(shù)，抽樣點數(shù)n1，補零后的數(shù)據(jù)長度，補零后的數(shù)據(jù)長度t2，抽，抽樣點數(shù)樣點數(shù)n2，則，則：0102112211 , ssffffntnt顯然，顯然，n2 n1，故故f2 f1。但并不代表：但并不代表：“”補零不能增加數(shù)據(jù)的有效長度，上面補零不能增加數(shù)據(jù)的有效長度，上面仍為仍為t1（有效抽樣點數(shù)（有效抽樣點數(shù)仍仍為為n1）。因此，。因此，補零是不能補零是不能提高頻率分辨率的。提高頻率分辨率的

10、。1.求出序列求出序列x(n)=cos(0.48n)+cos(0.52n)基于有限個樣點基于有限個樣點n=10的頻譜的頻譜;2.求求n=100時，取時，取x(n)的前的前10個，后個，后90個設(shè)為零，得到個設(shè)為零，得到x(n)的頻譜的頻譜;3.增加增加x(n)有效的樣點數(shù)，取有效的樣點數(shù)，取100個樣點得到個樣點得到x(n)的頻譜。的頻譜。n 基本思路基本思路n dit基基-2fftn 編程思想編程思想n 程序設(shè)計技術(shù)程序設(shè)計技術(shù)n 實序列的實序列的dft1010( )( )1( )( )nnknnnnknkx kx n wx nx k wnk=0, 1, , n-1 n=0, 1, , n-

11、1 復(fù)乘：復(fù)乘：2n復(fù)加：復(fù)加：) 1(nn實乘：實乘：24n實加：實加：2(21)nn 系數(shù)wnk具有以下特性： （1） 對稱性 nknnknww*)(（2） 周期性 )()(nknnknnnnknwww（3） 可約性 mnkmnnknnmkmnnknwwww/,（4）其它性質(zhì) ()()22/2(/2)2222cossin111n n kn n knknnnnjnnnk nknnjknknjknnjnnnnjnnnwwwwejwwweewee 2jnnwe1. 合并法：合并合并dft運算中的某些項。運算中的某些項。2. 分解法：將長序列將長序列dft利用對稱性和周期性，分解為短序利用對稱性和

12、周期性，分解為短序列列dft。 設(shè)序列x(n)長度為n，且滿足n=2m，m為正整數(shù)。按n的奇偶的奇偶把x(n)分解為兩個n/2點的子序列： 12, 1 , 0)() 12()()2(21nrrxrxrxrx若不滿足這個條件，可以若不滿足這個條件，可以人為地加上若干零值（加人為地加上若干零值（加零補長）使其達到零補長）使其達到 n=2mrknnrknrknnrkrnnrrknnrnnnnknnnnnnnknnknwrxwwrxwrxwrxwnxwnxwnxnxdftkx)()( )() 12()2()()()()()(2120221201)12(1202120101010為奇數(shù)為偶數(shù)則dft轉(zhuǎn)化

13、為下式：由于2/22222nnjnjnweew故上式可表示成：)()()()()(212/21202/1120kxwkxwrxwwrxkxknrknnrknrknnrrknnkrnww2/22/得到：)()()(212/120122/12011kxwrxwrxknxrknnrknrnnr同理： )(222kxknx說明：后半部分后半部分k值（值（n/2kn-1）所對應(yīng)的）所對應(yīng)的x1(k)、2(k)分別分別等于前半部分等于前半部分k值（值（0kn/2-1）所對應(yīng)的）所對應(yīng)的x1(k)、x2(k)。 即 我們已知： knknnnknnwwww2/2因此得：1221212( )( )( )( )(

14、 )222knnkknnx kx kw xknnnx kxkwxkx kw xk= -1x2(k)x1(k)knw1x1(k)knwx2(k)x1(k)knwx2(k)蝶形運算流圖按時間抽取將一個n點dft分解為兩個n/2點dft（n=8） x1(0)x1(1)x1(0) x(0)x(0)x(1)x1(2)x1(3)110nwdft2點nx1(1) x(2)x1(2) x(4)x1(3) x(6)x(2)x(3)x2(0)x2(1)x2(0) x(1)x(4)x(5)x2(2)x2(3)dft2點nx2(1) x(3)x2(2) x(5)x2(3) x(7)x(6)x(7)1nw2nw3nw1

15、1一個n點dft分解為兩個n/2點dft，每一個n/2點dft只需(n/2)2=n2/4次復(fù)數(shù)乘法，n/2(n/2-1)次復(fù)數(shù)加法。兩個n/2點dft共需2(n/2)2=n2/2次復(fù)數(shù)乘法和n(n/2-1)次復(fù)數(shù)加法。此外，把兩個n/2點dft合成為n點dft時，有n/2個蝶形運算，還需要n/2次復(fù)數(shù)乘法及2n/2=n次復(fù)數(shù)加法。因而通過第一步分解后，總共需要(n2/2)+(n/2)=n(n+1)/2n2/2次復(fù)數(shù)乘法和n(n/2-1)+n=n2/2次復(fù)數(shù)加法。由此可見，通過這樣分解后運算工作量差不多節(jié)省了一半。 既然這樣分解是有效的，由于n=2m，因而n/2仍是偶數(shù)，可以進一步把每個n/2點

16、子序列再按其奇偶部分分解為兩個n/4點的子序列。 依次類推，得到下圖所示的dft圖：dft4點nx3(0)x3(1)x3(0)x1(0)x(0)x3(1)x1(2)x(4)x1(0)x1(1)dft4點nx4(0)x4(1)x4(0)x1(1)x(2)x4(1)x1(3)x(6)x1(2)x1(3) 1 10nw2nwx(0)x(1)x(2)x(3)dft4點nx5(0)x5(1)x5(0)x2(0)x(1)x5(1)x2(2)x(5)x2(0)x2(1)dft4點nx6(0)x6(1)x6(0)x2(1)x(3)x6(1)x2(3)x(7)x2(2)x2(3) 1 1x(4)x(5)x(6)

17、x(7)0nw1nw2nw3nw 1 1 1 10nw2nwx(0)x(0)x(4)x(1) 10nwx(2)x(2)x(6)x(3) 10nw0nw2nw 1 1x(1)x(4)x(5)x(5) 10nwx(3)x(6)x(7)x(7) 10nw0nw2nw 1 10nw1nw2nw3nw 1 1 1 1由按時間抽取法fft的流圖可見，當(dāng)n=2m時，共有m級蝶形，每級都由n/2個蝶形運算組成，每個蝶形需要一次復(fù)乘、二次復(fù)加，因而每級運算都需n/2次復(fù)乘和n次復(fù)加，這樣m級運算總共需要：復(fù)乘數(shù) bnnnmabnnmnmff1122復(fù)加數(shù) 式中，數(shù)學(xué)符號：lb=log2。 由于計算機上乘法運算所

18、需的時間比加法運算所需的時間多得多，故以乘法為例，直接dft復(fù)數(shù)乘法次數(shù)是n2，fft復(fù)數(shù)乘法次數(shù)是(n/2) lbn。 直接計算dft與fft算法的計算量之比為 bnnbnnnmnn1212222當(dāng)n=2048時，這一比值為372.4，即直接計算dft的運算量是fft運算量的372.4倍。當(dāng)點數(shù)n越大時，fft的優(yōu)點更為明顯。 如果某通用單片計算機的速度為平均每次復(fù)數(shù)乘需要4s，每次復(fù)數(shù)加需要1s，用來計算n=1024點dft，問直接計算需要多少時間。用fft計算呢？照這樣計算，用fft進行快速卷積對信號進行處理時，估計可實現(xiàn)實時處理的信號最高頻率?！窘饨狻慨?dāng)n=1024=210時，直接計算

19、dft的復(fù)數(shù)乘法運算次數(shù)為：n 2 =10241024=1 048 576次復(fù)數(shù)加法運算次數(shù)為：n（n1）=10241023=1 047 552次所以，直接計算所用計算時間td為：td=1 048 576+1 047 552=5.241 856 s用fft計算1024點dft所需計算時間tf為：66f664 10lblb1 10210244 1010 1024 10 1 10230.72 msntnn n 快速卷積時，需要計算一次n點fft、n次頻域復(fù)數(shù)乘法和一次n點ifft。 所以，計算1024點快速卷積的計算時間tc約為：cf2102461440 s4 1024 s65536 stt次復(fù)數(shù)

20、乘計算時間所以， 每秒鐘處理的采樣點數(shù)（即采樣速率）s6102415 625 /65536 10f次 秒由采樣定理知， 可實時處理的信號最高頻率為smax156257.8125 khz22ff研讀教材研讀教材 p115的的4.2.4節(jié)。節(jié)。可以看出這種運算是很有規(guī)律的，其每級（每列）計算都是由n/2 個蝶形運算構(gòu)成的，每一個蝶形結(jié)構(gòu)完成下述基本迭代運算： rnmmmrnmmmwjxkxjxwjxkxkx)()()()()()(1111式中，m表示第m列迭代，k, j為數(shù)據(jù)所在行數(shù)。蝶形運算單元 1rnwxm 1(k)xm 1( j)xm(k) xm 1(k) xm 1( j)xm( j) xm

21、 1(k) xm 1( j)rnwrnw某一列的任何兩個節(jié)點k和j的節(jié)點變量進行蝶形運算后，得到結(jié)果為下一列k, j兩節(jié)點的節(jié)點變量，而和其他節(jié)點變量無關(guān)。因而可以采用原位運算，即某一列的n個數(shù)據(jù)送到存儲器后，經(jīng)蝶形運算，其結(jié)果為下一列數(shù)據(jù)，它們以蝶形為單位仍存以蝶形為單位仍存儲在這同一組存儲器中儲在這同一組存儲器中，直到最后輸出，中間無需其他存儲器。也就是蝶形的兩個輸出值仍放回蝶形的兩個輸入所在的存儲器中。每列的n/2個蝶形運算全部完成后，再開始下一列的蝶形運算。 這樣存儲器數(shù)據(jù)只需n個存儲單元。下一級的運算仍采用這種原位方式，只不過進入蝶形結(jié)的組合關(guān)系有所不同。這種原位運算結(jié)構(gòu)可以節(jié)省存儲

22、單元， 降低設(shè)備成本。 觀察同址計算結(jié)構(gòu)，發(fā)現(xiàn)當(dāng)運算完成后，fft的輸出x(k)按正常順序排列在存儲單元中，即按x(0)，x(1)，x(7)的順序排列，但是這時輸入x(n)卻不是按自然順序存儲的卻不是按自然順序存儲的，而是按x(0)，x(4)， ， x(7)的順序存入存儲單元，看起來好像是“混亂無序”的，實際上是有規(guī)律的，我們稱之為稱之為倒位序倒位序。 造成倒位序的原因是原因是輸入輸入x(n)按標(biāo)號按標(biāo)號n的偶奇的不斷分組的偶奇的不斷分組。 如果n用二進制數(shù)表示為(n2n1n0)2（當(dāng)n=8=23時，二進制為三位）， 第一次分組，n為偶數(shù)（相當(dāng)于n的二進制數(shù)的最低位n0=0）在上半部分，n為奇

23、數(shù)（相當(dāng)于n的二進制數(shù)的最低位 n0=1）在下半部分。下一次則根據(jù)次最低位n1為“0”或是“1”來分偶奇（而不管原來的子序列是偶序列還是奇序列而不管原來的子序列是偶序列還是奇序列）。如此繼續(xù)分下去，直到最后n個長度為1的子序列。下圖的樹狀圖描述了這種分成偶數(shù)子序列和奇數(shù)子序列的過程。 x(000)x(100)x(010)x(110)010101x(001)x(101)x(011)x(111)01010101x(n2n1n0)n2n1n0 一般實際運算中，總是先按自然順序?qū)⑤斎胄蛄写嫒氪鎯卧瑸榱说玫降刮恍虻呐帕?，我們通過變址運算來完成。如果輸入序列的自然順序號i用二進制數(shù)（例如n2n1n0）

24、表示，則其倒位序j對應(yīng)的二進制數(shù)就是（n0n1n2）。這樣，在原來自然順序時應(yīng)該放x(i)的單元，現(xiàn)在倒位序后應(yīng)放x(j)。 例如, n=8時，x(3)的標(biāo)號是i=3，它的二進制數(shù)是011，倒位序的二進制數(shù)是110，即j=6，所以原來存放在x(011)單元的數(shù)據(jù)現(xiàn)在應(yīng)該存放在x(110)內(nèi)。下表列出了n=8時的自然順序二進制數(shù)以及相應(yīng)的倒位序二進制數(shù)。 n=8時的自然順序二進制數(shù)和相應(yīng)的倒位序二進制數(shù)時的自然順序二進制數(shù)和相應(yīng)的倒位序二進制數(shù) 自然順序（自然順序（i） 二進制數(shù)二進制數(shù) 倒位序二進制數(shù)倒位序二進制數(shù) 倒位序（倒位序（j） 0123456700000101001110010111

25、011100010001011000110101111104261537 由表可見，自然順序數(shù)i增加1，是在順序數(shù)的二進制數(shù)最最低位加低位加1，向左進位向左進位。而倒序數(shù)j則是在二進制數(shù)最高位加最高位加1， 逢2向右進位向右進位。 例如，在（000）最高位加1，則得（100），再在（100）最高位加1，向右進位，則得（010）。因（100）最高位為1， 所以最高位加1要向次高位進位， 其實質(zhì)是將最高位變?yōu)?，再在次高位加1。用這種算法，可以從當(dāng)前任一倒序值求得下一個倒序值。 對于n=2m，m為二進制數(shù)最高位，其權(quán)值為n/2，且從左向右二進制位的權(quán)值依次為n/4，n/8， 2，1。 因此，最高位

26、加1相當(dāng)于十進制運算j+n/2。如果最高位是0（jn/2）， 則直接由j+n/2得下一個倒序值；如果最高位是1（jn/2），則要將最高位變?yōu)?（j j-n/2），次高位加1（j+n/4）。但次高位加1時，同樣要判斷次高位0、1值，如果為0（ji時，才將原存放x(i)及存放x(j)的存儲單元內(nèi)的內(nèi)容互換。這樣就得到輸入所需的倒位序列的順序。存儲單元自然順序輸入變址倒位序a(1)a(2)a(3)a(4)a(5)a(6)a(7)a(8)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)輸入是倒位序的，輸出是自然順序的。其第

27、一級（第一列）每個蝶形的兩節(jié)點間“距離”為1， 第二級每個蝶形的兩節(jié)點“距離”為2，第三級每個蝶形的兩節(jié)點“距離”為4。由此類推得，對n=2m點fft，當(dāng)輸入為倒位序，輸出為正常順序時，其第m級運算，每個蝶形的兩節(jié)點“距離”為2m-1。 由于對第m級運算，一個dft蝶形運算的兩節(jié)點“距離”為2m-1， 因而: rnmmmmmrnmmmmwkxkxkxwkxkxkx)2()()2()2()()(1111111 為了完成上兩式運算，還必須知道系數(shù)nr的變換規(guī)律： 把式中蝶形運算兩節(jié)點中的第一個節(jié)點標(biāo)號值，即k值，表示成m位（注意n=2m）二進制數(shù)； 把此二進制數(shù)乘上2m-m，即將此m位二進制數(shù)左移

28、m-m位（注意m是第m級運算），把右邊空出的位置補零， 此數(shù)即為所求r的二進制數(shù)。 從圖看出，wnr因子最后一列有n/2種，順序為wn0， wn1， ， 其余可類推。 12nnw 讀入讀入 x(n)，m n=2m 倒序倒序 l=1，m b=2l-1 j=0，b-1 p=2m-lj k=j，n-1，2l ( )()()( )()( )pnpnta ka kb wa kba ka kb wa kt 輸出輸出 function a=myfft_1(x,m) n=2m; if length(x)n x=x,zeros(1,n-length(x); end a=(x,n); for l=1:m b=2(l-1); nmr=2l; u=1; wn=exp(-i*2*pi/nmr); for j=1:b for k=j:nmr:n t=a(k+b)*u;