1229.FFT算法程序及分析

1、 數(shù)字信號(hào)處理 fft算法程序及分析摘 要：fft的算法程序分析主要分析了按時(shí)間抽?。╠it）的快速傅立葉變換的基2fft算法，通過(guò)對(duì)基2fft算法的原理的分析及與dft算法運(yùn)算量的比較，進(jìn)一步推導(dǎo)出了基rfft算法，重點(diǎn)是基rfft算法的推導(dǎo)。在具體的實(shí)例中，我們重點(diǎn)分析了fft過(guò)程中幅值大小與fft選用點(diǎn)數(shù)n的關(guān)系，驗(yàn)證fft變換的可靠性，考察在fft中數(shù)據(jù)樣本的長(zhǎng)度與dft的點(diǎn)數(shù)對(duì)頻譜圖的影響。關(guān)鍵字： 基2fft算法，基rfft算法，樣本長(zhǎng)度，選用點(diǎn)數(shù)要求：l 學(xué)習(xí)書(shū)上第六節(jié)的內(nèi)容，自己編程實(shí)現(xiàn)fft算法 。l 給出典型信號(hào)的時(shí)域和頻域圖，并加以分析。l 可嘗試實(shí)現(xiàn)分段卷積程序。l 論

2、文內(nèi)容含原程序、運(yùn)行結(jié)果，理論分析和典型信號(hào)時(shí)域圖。一快速傅立葉變換（fft）簡(jiǎn)介離散傅立葉變換（dft）是信號(hào)分析與處理中的一種重要的變換。因直接計(jì)算dft的計(jì)算量與變換區(qū)間長(zhǎng)度n的平方成正比，當(dāng)n較大時(shí)，計(jì)算量太大。所以在快速傅立葉變換（fft）出現(xiàn)以前，直接用dft算法進(jìn)行頻譜分析和信號(hào)的實(shí)時(shí)處理是不切實(shí)際的。1965年，庫(kù)利（j.w.cooley）和圖基（j.w.tukey）在計(jì)算數(shù)學(xué)雜志上發(fā)表了“機(jī)器計(jì)算傅立葉級(jí)數(shù)的一種算法”的文章，這是一篇關(guān)于計(jì)算dft的一種快速有效的計(jì)算方法的文章。它的思路建立在對(duì)dft運(yùn)算內(nèi)在規(guī)律的認(rèn)識(shí)之上。這篇文章的發(fā)表使dft的計(jì)算量大大減少，并導(dǎo)致了許多

3、計(jì)算方法的發(fā)現(xiàn)。這些算法統(tǒng)稱(chēng)為快速傅立葉變換（fast fourier transform），簡(jiǎn)稱(chēng)fft。1984年，法國(guó)的杜哈梅爾（p.dohamel）和霍爾曼（h.hollmann）提出的分裂基快速算法，使運(yùn)算效率進(jìn)一步提高?？焖俑盗⑷~變換（fft）不是一種新的變換，而是離散傅立葉變換（dft）的一種快速算法。fft分成兩大類(lèi)，即按時(shí)間抽?。╠ecimation-in-time，縮寫(xiě)為dit）法和按頻率抽?。╠ecimation-in-frequency，縮寫(xiě)為dif）法。fft算法主要包括基2fft算法，基4fft算法，混合基fft，基rfft算法和分裂基fft算法。二快速傅立葉變換（f

4、ft）定義 設(shè)x(n)為n點(diǎn)有限長(zhǎng)序列，其dft為 k=0,1,n-1 （1）其反變換（idft）為 n=0,1,n-1 （2）二者的差別只在于wn的指數(shù)符號(hào)不同，以及差一個(gè)常數(shù)因乘子1/n。一般x(n)和為復(fù)數(shù)序列。對(duì)某一個(gè)k值，直接按（1）式計(jì)算x(k)值需要n次復(fù)數(shù)乘法、（n-1）次復(fù)數(shù)加法。因此，對(duì)所有n個(gè)k值，共需要n2次復(fù)數(shù)乘法及n（n-1）次復(fù)數(shù)加法運(yùn)算。當(dāng)n1時(shí)，n（n-1）n2。（1）式可寫(xiě)為 所以，一次復(fù)數(shù)乘法需要四次實(shí)數(shù)乘法和二次實(shí)數(shù)加法；一次復(fù)數(shù)加法則需二次實(shí)數(shù)加法。因而每運(yùn)算一個(gè)x（k）需4n次實(shí)數(shù)乘法及2n+2（n-1）=2（2n-1）次實(shí)數(shù)加法。所以整個(gè)dft運(yùn)算

5、共需要4n2次實(shí)數(shù)乘法和n*2（2n-1）=2n（2n-1）次實(shí)數(shù)加法。由上述可見(jiàn)，n點(diǎn)dft的乘法和加法運(yùn)算次數(shù)均與n2成正比。當(dāng)n較大時(shí)，運(yùn)算量相當(dāng)可觀。所以，必須減少其運(yùn)算量，才能使dft在各種科學(xué)和工程計(jì)算中得到運(yùn)用。而dft運(yùn)算時(shí)間能否減少，關(guān)鍵在于時(shí)間運(yùn)算是否存在規(guī)律以及如何利用這些規(guī)律。仔細(xì)觀察dft的運(yùn)算可以看出，利用系數(shù)的一些固有特性，就可以減少運(yùn)算量：1) 的對(duì)稱(chēng)性 =2) 的周期性 =3） 的可約性 = =即有：= =-1 = 因此，利用這些性質(zhì)，可以合并dft運(yùn)算中的有些項(xiàng)；利用這些性質(zhì)可以將長(zhǎng)序列的dft分解為短序列的dft。從而，減少運(yùn)算次數(shù)，真正做到提高運(yùn)算效率。

6、三.fft基本形式1. 基2fft算法基2fft算法可以分為按時(shí)間抽?。╠it）的基-2fft算法（庫(kù)利-圖基算法）和按頻率抽取的基4fft算法。我們具體分析按時(shí)間抽?。╠it）的基-2fft算法。算法原理：先設(shè)序列點(diǎn)數(shù)為n=2l，l為整數(shù)。將n=2l的序列x(n)（n=0，1，n-1）先按n的奇偶分成兩組：， r=0,1,則可以將一個(gè)n點(diǎn)dft分解成兩個(gè)n/2點(diǎn)的dft，而x1 (r) 和x2(r)以及x1(k)和x2(k)都是n/2點(diǎn)的序列。x（k）表達(dá)為前后兩部分：前半部分 k=0,1, -1后半部分 k=0,1, -1只要求出0到（n/2-1）區(qū)間的所有x1(k)和x2(k)值，即可求

7、出0到（n-1）區(qū)間內(nèi)的所有x（k）的值。此時(shí)，我們可以利用蝶形信號(hào)流圖符號(hào)表示dft的運(yùn)算。下圖表示n=23=8的情況：n/2點(diǎn)dft n/2點(diǎn)dft -1 -1 -1 -1因?yàn)閚=2l，仍能被2整除。將x（k）分解后的兩部分按奇偶各分解為兩個(gè)點(diǎn)序列，而這兩個(gè)序列又可再分解為兩個(gè)點(diǎn)序列，依次類(lèi)推，可以一直分解到只有兩點(diǎn)的序列。由于n=2l,分解共需要l次。dft與fft運(yùn)算量的比較：由按時(shí)間抽選法fft的流圖可見(jiàn)，當(dāng)n=2l時(shí)，共有l(wèi)級(jí)蝶形，每級(jí)都由n/2個(gè)蝶形運(yùn)算組成，每個(gè)蝶形有一次復(fù)乘、二次復(fù)加，因此每級(jí)運(yùn)算都需要n/2次復(fù)乘和n次復(fù)加，這樣l級(jí)運(yùn)算總共需要復(fù)乘數(shù) 復(fù)加數(shù) 以乘法為例，直

8、接dft復(fù)數(shù)乘法次數(shù)是n2，fft復(fù)數(shù)乘法次數(shù)是。所以二者的計(jì)算量之比為2. 基r fft算法設(shè)序列x(n)長(zhǎng)度為n=rm, r和m均為任意整數(shù)。仿照基2dit-fft算法的推導(dǎo)方法，用m位r進(jìn)制數(shù)表示k和n 那么x（n）的dft可寫(xiě)成如下形式： （1-1）式中 為導(dǎo)出基r dit-fft算法，將n分解得 = 將上式代入（11）式得 （1-2）由上式同樣可得到m個(gè)遞推公式 （1-3）下面以n=42為例說(shuō)明基r dit-fft算法。將m=2,n4代入（1-3）式得基4 dit-fft遞推公式 (1-4)式中 （k1k0）和 (n1n0)分別表示k和n的m(m=2)位4進(jìn)制數(shù)。 表示4個(gè)4點(diǎn)dft

9、。 也表示4個(gè)4點(diǎn)dft。第一級(jí)的4點(diǎn)dft與第二級(jí)的4點(diǎn)dft通過(guò)旋轉(zhuǎn)因子連接起來(lái)，這與前面討論的基2dit-fft是一樣的。（1-4）式也可以寫(xiě)成矩陣形式。 x1(n0k0)的矩陣形式： =q式中q= 的矩陣形式： 完成4點(diǎn)dft的矩陣乘qa b c dt運(yùn)算的蝶形結(jié)符號(hào)如圖1（a）所示，對(duì)應(yīng)實(shí)際運(yùn)算流圖如圖1（b）所示。圖（b）只需要8次復(fù)加，但需要一次倒序。由x1(n0k0)和x2(k1k0)的矩陣表示式容易畫(huà)出16點(diǎn)基4dit-fft運(yùn)算流圖如圖2所示。 (a) 圖1（a）的蝶形結(jié)符號(hào) (b) 圖1（b）的實(shí)際運(yùn)算流圖 圖2 n16點(diǎn) 基4 dit-fft運(yùn)算流圖 從遞推公式（1-3

10、）的矩陣形式和圖都可看出，其輸入序列按二位四進(jìn)制倒序排列，而輸出x(k)則為順序排列。根據(jù)遞推公式（1-4）可寫(xiě)出基dit-fft算法的蝶形迭代公式 （1-5）式中 ，即 al(lrl+s)表示基r dit-fft算法第l級(jí)迭代中第l個(gè)運(yùn)算組的第s個(gè)輸出結(jié)點(diǎn)的中間結(jié)果。式（1-5）表明蝶形運(yùn)算流圖的輸入是按m位r進(jìn)制倒序，而輸出為順序。當(dāng)r=4,n=4m時(shí)，式（15）中 (1-6)式中，。推導(dǎo)上式的過(guò)程中用到關(guān)系式。將（1-6）式代入（1-5）式并按和展開(kāi)即可得到基算法的蝶形迭代公式 + + (1-7)式中 q=0,1,2, ,當(dāng)m=2,n=16時(shí),按（1-7）式畫(huà)出的流圖與圖2相同。如果對(duì)中

11、的k進(jìn)行分解，可導(dǎo)出基r dit-fft算法。限于篇幅，對(duì)此不進(jìn)行討論?；鵵 fft算法中，基4fft算法效率最高，觀察圖2可知，n點(diǎn) 基4 fft運(yùn)算流圖由m級(jí)運(yùn)算構(gòu)成，若不計(jì)乘j的運(yùn)算，則m級(jí)運(yùn)算中，從第二級(jí)開(kāi)始，每級(jí)要進(jìn)行n個(gè)旋轉(zhuǎn)因子的乘法運(yùn)算，所以復(fù)乘次數(shù)為 (1-8)其中包含乘的運(yùn)算。四.典型信號(hào)的分析fft的應(yīng)用1：已知信號(hào)由15hz幅值0.5的正弦信號(hào)和40hz幅值為2的正弦信號(hào)組成，數(shù)據(jù)采樣頻率為100hz，繪制n=128點(diǎn)的dft的幅頻圖和n=1024點(diǎn)的dft幅頻圖。即點(diǎn) 評(píng)：本題考察的是fft變換中，頻譜的對(duì)稱(chēng)性及幅值大小與fft選用點(diǎn)數(shù)n的關(guān)系。 源程序：clffs=1

12、00;n=128;n=0:n-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,n);mag=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(221)plot(f,mag);xlabel(frequency (hz);ylabel(magnitude);title(n=128)gridsubplot(222)plot(f(1:n/2),mag(1:n/2);xlabel(frequency (hz);ylabel(magnitude);title(n=128)gridfs=100;n=102

13、4;n=0:n-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,n);mag=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(223)plot(f,mag);xlabel(frequency (hz);ylabel(magnitude);title(n=1024)gridsubplot(224)plot(f(1:n/2),mag(1:n/2);xlabel(frequency (hz);ylabel(magnitude);title(n=128)grid程序運(yùn)行結(jié)果：總結(jié)：fs=100hz

14、，所以nyquist頻率為fs/2=50hz。由上圖可以看出：(a)、(c)整個(gè)頻譜圖都關(guān)于nyquist頻率對(duì)稱(chēng)。因此利用fft對(duì)信號(hào)做頻譜分析，只要考察0nyquist頻率（采樣頻率的一半）范圍的幅頻特性。再比較(a)(c)或(b)(d)可見(jiàn)，幅值大小和fft選用的點(diǎn)數(shù)n有關(guān)，但只要n足夠就可以不影響結(jié)果。fft的應(yīng)用2：用fft對(duì)信號(hào)進(jìn)行dft，對(duì)結(jié)果進(jìn)行ifft，比較idft的結(jié)果和原信號(hào)。點(diǎn) 評(píng)：本題旨在驗(yàn)證fft的可靠性。源程序： clf fs=100; %length of fft n=128; n=0:n-1; t=n/fs; x=sin(2*pi*40*t)+sin(2*pi

15、*15*t); subplot(221) plot(t,x); xlabel(t(s); ylabel(x); title(origianl signal) grid y=fft(x,n); mag=abs(y); f=(0:length(y)-1)*fs/length(y); subplot(222) plot(f(1:n/2),mag(1:n/2); xlabel(frequency (hz); ylabel(magnitude) title(fft to original signal) grid xifft=ifft(y); magx=real(xifft); ti=0:length(

16、xifft)-1/fs; subplot(223) plot(ti,magx); xlabel(frequency (hz); ylabel(magnitude) title(signal from ifft) grid yif=fft(xifft,n); mag=abs(yif); f=(0:length(y)-1)*fs/length(y); subplot(224) plot(f(1:n/2),mag(1:n/2); xlabel(frequency (hz); ylabel(magnitude) title(fft to signal from ifft) grid程序運(yùn)行結(jié)果：總結(jié)：

17、由圖可見(jiàn)，原信號(hào)和由ifft恢復(fù)信號(hào)的時(shí)域和頻域完全相同，因此fft可以由ifft恢復(fù)原信號(hào)。fft的應(yīng)用3：繪制信號(hào)(1)ndata=32,nfft=32;(2) ndata=32,nfft=128; (3) ndata=136,nfft=128; (4) ndata=136,nfft=512.的頻譜圖。點(diǎn) 評(píng)：本例主要是考察再fft中數(shù)據(jù)樣本的長(zhǎng)度與dft的點(diǎn)數(shù)對(duì)頻譜圖的影響。源程序： clf fs=100; %length of data ndata=32; %length of fft n=32; n=0:ndata-1; t=n/fs; x=0.5*sin(2*pi*15*t)+2*

18、sin(2*pi*40*t); y=fft(x,n); mag=abs(y); f=(0:length(y)-1)*fs/length(y); subplot(221) plot(f(1:n/2),mag(1:n/2); xlabel(frequency (hz); ylabel(magnitude) title(ndata=32 nfft=32) grid %length of data ndata=32; %length of fft n=128; n=0:ndata-1; t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); y=fft(x,n);

19、 mag=abs(y); f=(0:length(y)-1)*fs/length(y); subplot(222) plot(f(1:n/2),mag(1:n/2); xlabel(frequency (hz); ylabel(magnitude) title(ndata=32 nfft=128) grid %length of data ndata=136; %length of fft n=128; n=0:ndata-1; t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); y=fft(x,n); mag=abs(y); f=(0:length(

20、y)-1)*fs/length(y); subplot(223) plot(f(1:n/2),mag(1:n/2); xlabel(frequency (hz); ylabel(magnitude) title(ndata=136 nfft=128) grid %length of data ndata=136; %length of fft n=512; n=0:ndata-1; t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); y=fft(x,n); mag=abs(y); f=(0:length(y)-1)*fs/length(y); subp

21、lot(224) plot(f(1:n/2),mag(1:n/2); xlabel(frequency (hz); ylabel(magnitude) title(ndata=136 nfft=512) gri程序運(yùn)行結(jié)果：總結(jié)：對(duì)于y=ifft(x,n)，x為序列的傅立葉變換x(k) ;y是idftx(k),返回一個(gè)復(fù)數(shù)序列，；n為采用的idft的點(diǎn)數(shù)。當(dāng)n小于x長(zhǎng)度時(shí)，對(duì)x進(jìn)行截?cái)?；?dāng)n大于x的長(zhǎng)度時(shí)，對(duì)x進(jìn)行補(bǔ)零。因此，一般情況下，我們?nèi)〉狞c(diǎn)數(shù)和數(shù)據(jù)樣本相同，這樣頻譜具有較高的質(zhì)量，減小因補(bǔ)零或截?cái)喽a(chǎn)生的影響。五. 總 結(jié)通過(guò)以上的公式推導(dǎo)與例題分析，可以看出快速傅立葉變換的作用，fft并不是新的算法，是離散傅立葉變換的一種快速算法，正是因?yàn)橛辛怂攀筪ft的算法大大簡(jiǎn)化，在實(shí)際的信號(hào)分析處理中得到了廣泛的應(yīng)用。本次的科技小論文，僅僅是針對(duì)