快速傅里葉變換課件

1、數(shù)字信號處理AmplitudeTimeFrequency(a)數(shù)字信號處理AmplitudeTimeFrequency(a為了分析連續(xù)信號頻譜，需要使用DFT來逼近連續(xù)時間信號的傅里葉變換。但由于DFT要對連續(xù)信號進行采樣和截斷等處理，因此，常產(chǎn)生誤差分析，體現(xiàn)在三個方面。 混疊現(xiàn)象：由于不滿足抽樣定理而引起；截斷效應：由于序列截斷引起；柵欄效應：由于對頻譜的有限點采樣引起。信號逼近所產(chǎn)生的問題回顧2為了分析連續(xù)信號頻譜，需要使用DFT來逼近連續(xù)時間信號的傅回顧 頻譜混疊1. 混疊現(xiàn)象3回顧 頻譜混疊1. 混疊現(xiàn)象3回顧 頻譜混疊原因：時域采樣不滿足奈奎斯特準則： 其中 為抽樣頻率， 為信號最

2、高頻率。此條件只規(guī)定了 的下限，其上限要受頻率分辨率（抽樣間隔） 的約束。 也可表示為記錄長度的倒數(shù)，即：其中， 為抽樣點數(shù)。措施：采樣前對信號進行預濾波，采樣時滿足采樣定理，以避免信號在w=處附近的混迭，再濾去信號中頻率高于的頻率分量。4回顧 頻譜混疊原因：時域采樣不滿足奈奎斯特準則： 措施回顧 頻譜混疊需要注意的是：當 給定時, 若為了避免混疊現(xiàn)象而一味提高抽樣頻率 , 導致 增加，即導致頻率分辨能力下降; 反之, 若打算要提高頻率分辨率即減小 , 則導致減小 , 最終必須減小信號的高頻容量。顯然，針對高頻容量 與頻率分辨率 , 想保持其中一個不變而使另一個性能得以提高的唯一辦法是：增加記

3、錄長度內的點數(shù)N, 即：這是未采用任何特殊數(shù)據(jù)處理（例如加窗）情況下，為實現(xiàn)基本DFT算法所必須滿足的基本條件。5回顧 頻譜混疊需要注意的是：52. 截斷效應回顧 截斷效應窗函數(shù)不可能取無限寬，即其頻譜不可能為一沖激函數(shù)，信號的頻譜與窗函數(shù)的卷積必然產(chǎn)生拖尾、變寬的現(xiàn)象，從而造成截斷效應 。如下圖所示。62. 截斷效應回顧 截斷效應窗函數(shù)不可能取無限寬，即“泄漏”是由矩形窗函數(shù)帶來的，求極限：的連續(xù)頻譜。為處的一根譜線，變成了以原來在為中心，形狀可見：回顧 截斷效應7“泄漏”是由矩形窗函數(shù)帶來的，求極限：的連續(xù)頻譜。為處的一根小、主瓣窄的窗函數(shù)。接近為了盡量減少泄漏，需要尋找頻域中窗函數(shù)具有旁

4、瓣，并且主瓣也不過，泄漏的產(chǎn)生是由于所以不能用無限加寬窗口來減少泄漏。，意味著需要無限加寬窗函數(shù)，等于未截短。若使泄漏至0，則回顧 截斷效應占有一定寬度。，即旁瓣8小、主瓣窄的窗函數(shù)。接近為了盡量減少泄漏，需要尋找頻域中窗函回顧 柵欄效應N點DFT是在區(qū)間0, 2上的N點等間隔采樣，采樣點之間的頻譜函數(shù)值是不知道的，就好像從N+1個柵欄縫隙中觀看信號的頻譜特性，得到的是N個縫隙中看到的頻譜函數(shù)值，這種現(xiàn)象稱為柵欄效應。原因：對信號的頻譜進行有限點采樣；后果：柵欄效應可能漏掉（擋?。┐蟮念l譜分量；措施：對原序列補0，增大N，以增加采樣點。3. 柵欄效應管中窺豹9回顧 柵欄效應N點DFT是在區(qū)間0

5、, 2上的N減小柵欄效應的一個方法是，在所取數(shù)據(jù)的末端添加一些零值點，使一個周期內點數(shù)增加, 但是不改變原有的記錄數(shù)據(jù)。這種方法相當于加長了周期T0，T0增加, 使得頻域抽樣間隔變小, 從而能保持原來頻譜形式不變的情況下使譜線變密, 也就使頻譜抽樣點數(shù)增加。這樣，原來看不到的頻譜分量就有可能看到了?；仡?柵欄效應 減小柵欄效應的方法 補零10減小柵欄效應的一個方法是，在所取數(shù)據(jù)的末端添加一些零值點，回顧 柵欄效應8176542308111654230109712后面補4個零值，則有：11回顧 柵欄效應81765423081116542301【解】例 1： 柵欄效應已知：12【解】例 1： 柵欄

6、效應已知：12例 1： 柵欄效應x=2,3,3,2;subplot(2,2,1);stem(x,fill);title(x(n);N=8*length(x);NFFT = 2nextpow2(N);Xk=fft(x,NFFT);subplot(2,2,2);f = 0:2*pi/(NFFT-1):2*pi;stem(f,abs(Xk);title(X(k) DFT點數(shù) N= num2str(NFFT);x1=2,3,3,2,0,0,0,0;subplot(2,2,3);stem(x1,fill);title(x1(n);N1=8*length(x1);NFFT1 = 2nextpow2(N1)

7、;X1k=fft(x1,NFFT1);subplot(2,2,4);f1 = 0:2*pi/(NFFT1-1):2*pi;stem(f1,abs(X1k);title(X1(k) DFT點數(shù) N= num2str(NFFT1);13例 1： 柵欄效應x=2,3,3,2;13注意：補加零點以改變周期時, 所用窗函數(shù)寬度卻不能變。亦即必須按數(shù)據(jù)記錄原長來選擇窗函數(shù), 而不能按補了零值點后的長度來選擇窗函數(shù)。即：“先加窗, 再補零。”回顧 柵欄效應14注意：回顧 柵欄效應14一般來說，信號長度 越長，抽樣點數(shù) 越大，則 越小，即分辨力越強。但此處的 是指實際的信號長度， 是指這個長度上的抽樣點數(shù)，而

8、不是進行補零等輔助措施后的信號長度及抽樣點數(shù)。回顧 分辨率 頻率分辨率所謂頻率分辨率，是指能夠將信號中兩個靠得很近的譜峰分開的能力：15一般來說，信號長度 越長，抽樣點數(shù) 越由上兩式可見，當采樣頻率由變?yōu)闀r，在數(shù)據(jù)長度變?yōu)?，頻率采樣間隔不變。不變情況下，采樣點數(shù)由僅增加采樣率并不能改變頻率分辨率。如果數(shù)據(jù)長度不變，則有： 注意（1）：單純增加采樣率并不能改變頻率分辨率 ！回顧 分辨率分析：顯然，只有“增加點數(shù)的同時導致數(shù)據(jù)有效長度的增加”才能使分辨率越好。 16由上兩式可見，當采樣頻率由變?yōu)闀r，在數(shù)據(jù)長度變?yōu)?，頻率采樣間例 2： 分辨率已知某信號如下：假設最高頻率 ，請用 對其抽樣，并觀察其頻

9、譜情況。【解】 假設取信號長度為： ，經(jīng)抽樣所得的 的點數(shù)為：對 做DFT變換，所得到的頻率最大分辨率為：17例 2： 分辨率已知某信號如下：【解】 17例 2： 分辨率clear all;L=256; % 信號長度N=256; % DFT點數(shù)n=0:1:L-1;f1=2;f2=2.02;f3=2.07x=sin(2*pi*f1*0.1*n)+sin(2*pi*f2*0.1*n)+sin(2*pi*f3*0.1*n);subplot(2,1,1);stem(x);title(x(n) 信號長度 tp= num2str(L);X=fft(x,N);Xk=abs(X(1:N/2);k=0:1:N/

10、2-1;w=2*pi/N*k;subplot(2,1,2);plot(5*w/pi,Xk/N);title(DFT 點數(shù)N= num2str(N);axis(1,3,0,1);現(xiàn)象： 只顯示了兩個頻率：f1、f3原因： f2-f1=0.02F0，所以能分辨出來改進方法： 增加采樣點數(shù) N，即增加數(shù)據(jù)長度。18例 2： 分辨率clear all;現(xiàn)象：18例 2： 分辨率clear all;L=1024; % 信號長度N=1024; % DFT點數(shù)19例 2： 分辨率clear all;19注意（2）：并非補零后就能提高頻率分辨率 ！設原數(shù)據(jù)長度T1，抽樣點數(shù)N1，補零后的數(shù)據(jù)長度T2，抽樣點數(shù)

11、N2，則：顯然，N2 N1，故F2 F1。但并不代表：“補零后，頻率分辨率提高了！”因為：補零不能增加數(shù)據(jù)的有效長度，上面實際數(shù)據(jù)的有效長度仍為T1（有效抽樣點數(shù)仍為N1）。因此，補零是不能提高頻率分辨率的。回顧 分辨率20注意（2）：并非補零后就能提高頻率分辨率 ！顯然，N2 例 3： 分辨率求出序列x(n)=cos(0.48n)+cos(0.52n)基于有限個樣點n=10的頻譜;求n=100時，取x(n)的前10個，后90個設為零，得到x(n)的頻譜;增加x(n)有效的樣點數(shù)，取100個樣點得到x(n)的頻譜。21例 3： 分辨率求出序列x(n)=cos(0.48n)例 3： 分辨率x(n

12、)基于10個樣點的頻譜figure(1)n=0:1:99;x=cos(0.48*pi*n)+cos(0.52*pi*n);n1=0:1:9;y1=x(1:1:10);subplot(2,1,1);stem(n1,y1);title(signal x(n), 0 = n = 9);xlabel(n)axis(0,10,-2.5,2.5)Y1=fft(y1);magY1=abs(Y1(1:1:6);k1=0:1:5;w1=2*pi/10*k1;subplot(2,1,2);stem(w1/pi,magY1);title(10點DFT);xlabel(w/pi)，axis(0,1,0,10) 22例

13、 3： 分辨率x(n)基于10個樣點的頻譜22例 3： 分辨率在10個樣點的基礎上添90個零，得到密度高的頻譜figure(2)n3=0:1:99;y3=x(1:1:10) zeros(1,90); 添90個零。得到100個數(shù)據(jù)subplot(2,1,1);stem(n3,y3);title(signal x(n), 0 = n = 9 + 90 zeros);xlabel(n)axis(0,100,-2.5,2.5)Y3=fft(y3);magY3=abs(Y3(1:1:51);k3=0:1:50;w3=2*pi/100*k3;subplot(2,1,2);stem(w3/pi,magY3)

14、;title(100點DFT);xlabel(w/pi)axis(0,1,0,10)23例 3： 分辨率在10個樣點的基礎上添90個零，得到密度例 3： 分辨率增加x(n)有效的樣點數(shù)，取100個樣點figure(3)n=0:1:99;x=cos(0.48*pi*n)+cos(0.52*pi*n);subplot(2,1,1);stem(n,x);title(signal x(n), 0 = n = 99);xlabel(n)axis(0,100,-2.5,2.5)X=fft(x);magX=abs(X(1:1:51);k=0:1:50;w=2*pi/100*k;subplot(2,1,2);

15、stem(w/pi,magX);title(100點DFT);xlabel(w/pi)axis(0,1,0,60) 24例 3： 分辨率增加x(n)有效的樣點數(shù)，取100個樣點4 快速傅里葉變換（FFT） 基本思路 DIT基-2FFT 編程思想 程序設計技術 實序列的DFT254 快速傅里葉變換（FFT） 基本思路25圖基 庫利FFT引言 “機器計算傅里葉級數(shù)的一種算法”傅里葉FT26圖基 一、問題的提出k=0, 1, , N-1 n=0, 1, , N-1 4.1 基本思路計算量復乘：復加：實乘：實加：DFT都是復數(shù)27一、問題的提出k=0, 1, , N-1 n=0, 1, 二、改善原理系

16、數(shù)Wnk具有以下特性： （1） 對稱性 （2） 周期性 4.1 基本思路（3） 可約性 （4）其它性質 28二、改善原理系數(shù)Wnk具有以下特性： （2） 周期性 4合并法：合并DFT運算中的某些項。分解法：將長序列DFT利用對稱性和周期性，分解為短序列DFT。三、改善DFT運算效率的基本途徑4.1 基本思路快速傅里葉變換正是基于這樣思想而發(fā)展起來的。它的算法形式有很多種，但基本上可以分成兩大類：按時間抽?。―ecimation in Time，DIT）按頻率抽取（Decimation in Frequency， DIF）29合并法：合并DFT運算中的某些項。三、改善DFT運算效率的基一、算法原

17、理 設序列x(n)長度為N，且滿足N=2M，M為正整數(shù)。按n的奇偶把x(n)分解為兩個N/2點的子序列： 4.2 按時間抽取的基 -2 FFT若不滿足這個條件，可以人為地加上若干零值（加零補長）使其達到 N=2M則DFT轉化為下式：30一、算法原理4.2 按時間抽取的基 -2 FFT若不滿足這個由于4.2 按時間抽取的基 -2 FFT故上式可表示成：可見，一個N點DFT分解成兩個N/2點的DFT。但由于X1(k)，X2(k)只有N/2個點，而X(k)卻有N個點，故計算得到的只是X(k)的前一半的結果，要用X1(k)、X2(k)來表達全部的X(k)值，還必須應用系數(shù)的周期性。31由于4.2 按時

18、間抽取的基 -2 FFT故上式可表示成：可見得到：同理： 說明：后半部分k值（N/2kN-1）所對應的X1(k)、2(k)分別等于前半部分k值（0kN/2-1）所對應的X1(k)、X2(k)。 4.2 按時間抽取的基 -2 FFT即 32得到：同理： 說明：后半部分k值（N/2kN-1）所對應我們已知： 因此得：4.2 按時間抽取的基 -2 FFT= -1蝶形運算流圖33我們已知： 因此得：4.2 按時間抽取的基 -2 FFT= 按時間抽取將一個N點DFT分解為兩個N/2點DFT（N=8） 4.2 按時間抽取的基 -2 FFT34按時間抽取將一個N點DFT分解為兩個N/2點DFT（N=8）一個

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

20、再按其奇偶部分分解為兩個N/4點的子序列。 依次類推，得到下圖所示的DFT圖：4.2 按時間抽取的基 -2 FFT35一個N點DFT分解為兩個N/2點DFT，每一個N/2點DF一個N=8點DFT分解為四個N/4點DFT 4.2 按時間抽取的基 -2 FFT36一個N=8點DFT分解為四個N/4點DFT 4.2 按時間抽N=8 按時間抽取的蝶形運算流圖4.2 按時間抽取的基 -2 FFT37N=8 按時間抽取的蝶形運算流圖4.2 按時間抽取的基 -2二、運算量比較由按時間抽取法FFT的流圖可見，當N=2M時，共有M級蝶形，每級都由N/2個蝶形運算組成，每個蝶形需要一次復乘、二次復加，因而每級運算

21、都需N/2次復乘和N次復加，這樣M級運算總共需要：復乘數(shù) 復加數(shù) 式中，數(shù)學符號：lb=log2。 4.2 按時間抽取的基 -2 FFT38二、運算量比較復乘數(shù) 復加數(shù) 式中，數(shù)學符號：lb=log2由于計算機上乘法運算所需的時間比加法運算所需的時間多得多，故以乘法為例，直接DFT復數(shù)乘法次數(shù)是N2，F(xiàn)FT復數(shù)乘法次數(shù)是(N/2) lbN。 直接計算DFT與FFT算法的計算量之比為 當N=2048時，這一比值為372.4，即直接計算DFT的運算量是FFT運算量的372.4倍。當點數(shù)N越大時，F(xiàn)FT的優(yōu)點更為明顯。 4.2 按時間抽取的基 -2 FFT39由于計算機上乘法運算所需的時間比加法運算

22、所需的時間多得多，故如果某通用單片計算機的速度為平均每次復數(shù)乘需要4s，每次復數(shù)加需要1s，用來計算N=1024點DFT，問直接計算需要多少時間。用FFT計算呢？照這樣計算，用FFT進行快速卷積對信號進行處理時，估計可實現(xiàn)實時處理的信號最高頻率。習題 4.1【解】當N=1024=210時，直接計算DFT的復數(shù)乘法運算次數(shù)為：N 2 =10241024=1 048 576次復數(shù)加法運算次數(shù)為：N（N1）=10241023=1 047 552次所以，直接計算所用計算時間TD為：TD=41061 048 576+1 047 5521106=5.241 856 s用FFT計算1024點DFT所需計算時

23、間TF為：40如果某通用單片計算機的速度為平均每次復數(shù)乘需要4s，每次快速卷積時，需要計算一次N點FFT、N次頻域復數(shù)乘法和一次N點IFFT。 所以，計算1024點快速卷積的計算時間Tc約為：習題 4.141快速卷積時，需要計算一次N點FFT、N次頻域復數(shù)乘法和一次所以， 每秒鐘處理的采樣點數(shù)（即采樣速率）由采樣定理知， 可實時處理的信號最高頻率為習題 4.142所以， 每秒鐘處理的采樣點數(shù)（即采樣速率）由采樣定理知， 研讀教材 P115的4.2.4節(jié)。一、原位運算（同址運算）可以看出這種運算是很有規(guī)律的，其每級（每列）計算都是由N/2 個蝶形運算構成的，每一個蝶形結構完成下述基本迭代運算：

24、式中，m表示第m列迭代，k, j為數(shù)據(jù)所在行數(shù)。4.3 編程思想43研讀教材 P115的4.2.4節(jié)。一、原位運算（同址運算）式蝶形運算單元 4.3 編程思想44蝶形運算單元 4.3 編程思想44某一列的任何兩個節(jié)點k和j的節(jié)點變量進行蝶形運算后，得到結果為下一列k, j兩節(jié)點的節(jié)點變量，而和其他節(jié)點變量無關。因而可以采用原位運算，即某一列的N個數(shù)據(jù)送到存儲器后，經(jīng)蝶形運算，其結果為下一列數(shù)據(jù)，它們以蝶形為單位仍存儲在這同一組存儲器中，直到最后輸出，中間無需其他存儲器。也就是蝶形的兩個輸出值仍放回蝶形的兩個輸入所在的存儲器中。每列的N/2個蝶形運算全部完成后，再開始下一列的蝶形運算。 這樣存儲

25、器數(shù)據(jù)只需N個存儲單元。下一級的運算仍采用這種原位方式，只不過進入蝶形結的組合關系有所不同。這種原位運算結構可以節(jié)省存儲單元， 降低設備成本。 4.3 編程思想45某一列的任何兩個節(jié)點k和j的節(jié)點變量進行蝶形運算后，得到結果二、倒位序規(guī)律 觀察同址計算結構，發(fā)現(xiàn)當運算完成后，F(xiàn)FT的輸出X(k)按正常順序排列在存儲單元中，即按X(0)，X(1)，X(7)的順序排列，但是這時輸入x(n)卻不是按自然順序存儲的，而是按x(0)，x(4)， ， x(7)的順序存入存儲單元，看起來好像是“混亂無序”的，實際上是有規(guī)律的，我們稱之為倒位序。 4.3 編程思想46二、倒位序規(guī)律4.3 編程思想46 造成倒

26、位序的原因是輸入x(n)按標號n的偶奇的不斷分組。 如果n用二進制數(shù)表示為(n2n1n0)2（當N=8=23時，二進制為三位）， 第一次分組，n為偶數(shù)（相當于n的二進制數(shù)的最低位n0=0）在上半部分，n為奇數(shù)（相當于n的二進制數(shù)的最低位 n0=1）在下半部分。下一次則根據(jù)次最低位n1為“0”或是“1”來分偶奇（而不管原來的子序列是偶序列還是奇序列）。如此繼續(xù)分下去，直到最后N個長度為1的子序列。下圖的樹狀圖描述了這種分成偶數(shù)子序列和奇數(shù)子序列的過程。 4.3 編程思想47 造成倒位序的原因是輸入x(n)按標號n的偶奇倒位序的形成 4.3 編程思想48倒位序的形成 4.3 編程思想48 一般實際

27、運算中，總是先按自然順序將輸入序列存入存儲單元，為了得到倒位序的排列，我們通過變址運算來完成。如果輸入序列的自然順序號I用二進制數(shù)（例如n2n1n0）表示，則其倒位序J對應的二進制數(shù)就是（n0n1n2）。這樣，在原來自然順序時應該放x(I)的單元，現(xiàn)在倒位序后應放x(J)。 例如, N=8時，x(3)的標號是I=3，它的二進制數(shù)是011，倒位序的二進制數(shù)是110，即J=6，所以原來存放在x(011)單元的數(shù)據(jù)現(xiàn)在應該存放在x(110)內。下表列出了N=8時的自然順序二進制數(shù)以及相應的倒位序二進制數(shù)。 4.3 編程思想49 一般實際運算中，總是先按自然順序將輸入序列N=8時的自然順序二進制數(shù)和相

28、應的倒位序二進制數(shù) 自然順序（I） 二進制數(shù) 倒位序二進制數(shù) 倒位序（J） 01234567000001010011100101110111000100010110001101011111042615374.3 編程思想50N=8時的自然順序二進制數(shù)和相應的倒位序二進制數(shù) 自然順序（ 由表可見，自然順序數(shù)I增加1，是在順序數(shù)的二進制數(shù)最低位加1，向左進位。而倒序數(shù)J則是在二進制數(shù)最高位加1， 逢2向右進位。 例如，在（000）最高位加1，則得（100），再在（100）最高位加1，向右進位，則得（010）。因（100）最高位為1， 所以最高位加1要向次高位進位， 其實質是將最高位變?yōu)?，再在次高

29、位加1。用這種算法，可以從當前任一倒序值求得下一個倒序值。 4.3 編程思想51 由表可見，自然順序數(shù)I增加1，是在順序數(shù)的二 對于N=2M，M為二進制數(shù)最高位，其權值為N/2，且從左向右二進制位的權值依次為N/4，N/8， 2，1。 因此，最高位加1相當于十進制運算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)的存儲單元內的內容互換。這樣就得到輸入所需的倒位序列的順序。4.3 編程

30、思想53 把按自然順序存放在存儲單元中的數(shù)據(jù)，換成FFN=8 倒位序的變址處理 4.3 編程思想54N=8 倒位序的變址處理 4.3 編程思想543. 蝶形運算兩節(jié)點的“距離”輸入是倒位序的，輸出是自然順序的。其第一級（第一列）每個蝶形的兩節(jié)點間“距離”為1， 第二級每個蝶形的兩節(jié)點“距離”為2，第三級每個蝶形的兩節(jié)點“距離”為4。由此類推得，對N=2M點FFT，當輸入為倒位序，輸出為正常順序時，其第m級運算，每個蝶形的兩節(jié)點“距離”為2m-1。 4.3 編程思想553. 蝶形運算兩節(jié)點的“距離”4.3 編程思想55四、WNr 的確定 由于對第m級運算，一個DFT蝶形運算的兩節(jié)點“距離”為2m

31、-1， 因而: 為了完成上兩式運算，還必須知道系數(shù)Nr的變換規(guī)律： 把式中蝶形運算兩節(jié)點中的第一個節(jié)點標號值，即k值，表示成M位（注意N=2M）二進制數(shù)； 4.3 編程思想56四、WNr 的確定 為了完成上兩式運算，還必須 把此二進制數(shù)乘上2M-m，即將此M位二進制數(shù)左移M-m位（注意m是第m級運算），把右邊空出的位置補零， 此數(shù)即為所求r的二進制數(shù)。 從圖看出，WNr因子最后一列有N/2種，順序為WN0， WN1， ， 其余可類推。 4.3 編程思想57 把此二進制數(shù)乘上2M-m，即將此M位二進4.4 FFT程序設計程序框圖 function A=myFFT_1(x,m) N=2m; if

32、length(x)N% 補0 x=x,zeros(1,N-length(x); end A=myBitRevorder(x,N); % 倒序函數(shù) for L=1:m % L級循環(huán) B=2(L-1); % B個旋轉因子 Nmr=2L; u=1; % 旋轉因子u初始化 WN=exp(-i*2*pi/Nmr); % 因子WN for j=1:B % B次蝶形運算 for k=j:Nmr:N% 本次碟形運算 t=A(k+B)*u;% 碟形運算 A(k+B)=A(k)-t; A(k)=A(k)+t; end u=u*WN; % 修改旋轉因子 end end584.4 FFT程序設計程序框圖 function A=myFfunction A=my