手把手教你理解(FFT)

手把手教你理解(FFT)實(shí)驗(yàn)原理傅里葉變換是一種將信號(hào)從時(shí)域變換到頻域的變換形式，是信號(hào)處理的重要分析工具。離散傅里葉變換(DFT)是傅里葉變換在離散系統(tǒng)中的表示形式。但是DFT的計(jì)算量非常大，F(xiàn)FT就是DFT的一種快速算法,FFT將DFT的N?步運(yùn)算減少至(N/2)log?N步。離散信號(hào)x(n)的傅里葉變換可以表示為：N-1X(R)=。［觀jV=O式中的Wn稱為蝶形因子，利用它的對(duì)稱性和周期性可以減少運(yùn)算量。一般而言,FFT算法分為時(shí)間抽取(DIT)和頻率抽取(DIF)兩大類。兩者的區(qū)別是蝶形因子出現(xiàn)的位置不同，前者中蝶形因子出現(xiàn)在輸入端，后者中出現(xiàn)在輸出端。本實(shí)驗(yàn)以時(shí)間抽取方法為例。前者和后者如圖所示：時(shí)間抽取FFT是將N點(diǎn)輸入序列x(n)按照偶數(shù)項(xiàng)和奇數(shù)項(xiàng)分解為偶序列和奇序列。偶序列為：x(0),x(2),x(4),...,x(N-2)；奇序列為：x(l),x(3),x(5),...,x(N-l)o這樣x(n)的N點(diǎn)DFT可寫(xiě)成：N/2-1 jV/2-1X（k）=￡.伽）叱嚴(yán)+￡x（2〃+l）W嚴(yán)於TOC\o"1-5"\h\z”=0 "=0考慮到Wn的性質(zhì)，即叱=［嚴(yán)”竹2=嚴(yán)/(N⑵二％因此有：N/2-1 N/2-1X伙)=￡x(2〃)w為+叱￡x(2〃+l)叱2H=0 /1=0或者寫(xiě)成：X伙)二瞅)+W剛由于Y(k)與Z(k)的周期為N/2,并且利用Wn的對(duì)稱性和周期性，即:W嚴(yán)4=_wk可得：X(k+N/2)二除)—W；Z(k)對(duì)Y(k)與Z(k)繼續(xù)以同樣的方式分解下去，就可以使一個(gè)N點(diǎn)的DFT最終用一組2點(diǎn)的DFT來(lái)計(jì)算。在基數(shù)為2的FFT中，總共有l(wèi)og2(N)級(jí)運(yùn)算，每級(jí)中有N/2個(gè)2點(diǎn)FFT蝶形運(yùn)算。單個(gè)蝶形運(yùn)算示意圖如下:Wn-1

Wn-1以N=8為例，時(shí)間抽取FFT的信號(hào)流圖如下:x(0)x(4)xx(0)x(4)x⑵x(6)x(l)x(5)x⑶x(7)X(0)X(l)X⑵X⑶X⑷X⑸X⑹X⑺從上圖可以看出，輸出序列是按自然順■列的,而輸入序列的順序則是“比特反轉(zhuǎn)”方式排列的。也就是說(shuō)，將字號(hào)用二進(jìn)制表示，然后將二進(jìn)制數(shù)以相反方向排列，再以這個(gè)數(shù)作為序號(hào)。如011變成110,那么第3個(gè)輸入值和第六個(gè)輸入值就要交換位置了。本實(shí)驗(yàn)中釆用了一種比較常用有效的方法完成這一步工作雷徳算法。

/*PiogiamforFFT*/#include<math.h>#include<stdio.h>#deflneN64constfloatpi=3.1415926;/*FFT點(diǎn)數(shù)*//*輸入信號(hào)序列/*輸入信號(hào)序列*//*輸出頻譜序列*/floatv_ie[64],y_im[64];floatvout[64];/*蝶形因子/*蝶形因子*//*蝶形運(yùn)算的級(jí)數(shù)，即Log2(N)*//*臨時(shí)變量*/iiitm;floatt_re,t_im,v_re?v_im;iiita,b,c;voidpaixu(floaty_re[N].floaty_im[N])/*用雷徳算法對(duì)輸入信號(hào)序列進(jìn)行倒序重排*/尸0；foi(i=0;i<N;i++)t_re=y_ie|j];y-ieU]=y_re[i];y_imU]=y_mi[i];y_re[i]=t_ie;}k=N/2;while((k<=^j)&(k>0)){j=j-k；k=k/2;

尸j+k;voidfft(floaty_ie[N].floatv_im[N]){/*計(jì)算蝶形運(yùn)算的級(jí)數(shù)log2(N)*/f=N;foi(m=l;(仁f/2)!=l;m十十);t_ie=v_ie*w_ie-v_mi*w_ini;t_ini=v_re*w_im+v_mi*w_ie;v_ie=t_ie;***fft***/***fft***/foi(n=l;n<=m;n++){a=l; /*a=2foi(n=l;n<=m;n++){a=l; /*a=2的n次方*/foi(i=0;i<n;i++)a=a*2;b=a/2;v_ie=1.0; /*蝶形因子*/v_un=0.0;w_re=cos(pi/b);w_im=-sin(pi/b);foi(j=0;j<b;j卄)/*蝶形運(yùn)算*/c=i十b;t_ie=v_ie[c]*v_ie-y_nn[c]t_im=y_ie[c]*v_im+y_im[c]*v_ie;y_re[c]=y_ie[i]-t_re;v_re[i]=y_ie[i]+t_re;v_un=t_im;}}}voidmaui(){/*初始化數(shù)據(jù)空間*/for(i=0;i<N;i十十){x_re[i]=0;x_im[i]=0;}〃輸入的模擬信號(hào)for(i=0;i<=N;i++){//x_ie[i]=(2*cos(-16*pi*i/N)+cos(-10*pi*i/N)+5*cos(-24*pi*i/N));x_re[i]=4*cos(-10*pi*i/N)+2*cos(-16*pi*i/N);x_un[i]=0;}/*復(fù)制到輸出數(shù)組*/for(i=0;i<N;i十十){y_ie[i]=x_ie[i];y_mi[i]=x_mi[i];}paixu(v_re,y_im); 〃對(duì)輸入數(shù)據(jù)進(jìn)行倒序fft(y_re,y_im); 〃對(duì)輸入數(shù)據(jù)的點(diǎn)數(shù)進(jìn)行FFT運(yùn)算for(i=0;i<N/2;i-H-){yout[i]=sqrt(y_ie[i]*y_re[i]+y_im[i]*v_mi[i])*2/N; 〃計(jì)算幅值while(l);}實(shí)驗(yàn)步驟以64點(diǎn)FFT的信號(hào)流圖為例，理解FFT算法的過(guò)程；在CCS3.3環(huán)境中打開(kāi)本實(shí)驗(yàn)的工程(Example_fflpjt),編譯并重建.out輸出文件，然后通過(guò)仿真器把執(zhí)行代碼下載到DSP芯片中；運(yùn)行程序；選擇view->grapli->time/fiequency...。設(shè)置對(duì)話框中的參數(shù)：其中“StartAddress設(shè)為"x_re”，"Acquisitionbuffersize"和“DisplayDatasize”都設(shè)為“64"，并且把“DSPDataType”設(shè)為“32-bit floatuigpoint”(如圖)，SSGraphPropertyDialogDisplayTypeSingleTimeAGraphTitleGraphicalDisplayStartAddressx_rePageDataAcq,uisitionBufferSize64言*Irkereirient1DisplayD：at色Size64DSPDataType32"bitfloatingpointSajnpliiigRateCHz)1PlotDataPromLefttoRightLefJisplayYesAutoscaleOnDCValue0AxesDisplOnTimeDisplayUnit5StatusBarDisplayOnV|QK|QancJ|Help|同樣方法觀察經(jīng)DFT變換后的輸出序列“y」e”的波形，“StartAddress"改為“y」e”，其余參數(shù)不變（如圖）;釆樣定理告訴我們，釆樣頻率要大于信號(hào)頻率的兩借，假設(shè)采樣頻率為Fs,信號(hào)頻率F,釆樣點(diǎn)數(shù)為N。模糊知識(shí)點(diǎn)的解釋：1?采樣出的點(diǎn)經(jīng)過(guò)FFT運(yùn)算后，各自點(diǎn)對(duì)應(yīng)的頻率是多少？答：也可以看做是將第一個(gè)點(diǎn)分做兩半分，另一半移到最后）則表示釆樣頻率Fs,這中間被N-1個(gè)點(diǎn)平均分成N等份，每個(gè)點(diǎn)的頻率依次增加。例如某點(diǎn)n所表示的頻率為：Fn=（n-l）*Fs/No舉例：Fn所能分辨到頻率為為Fs/N,如果釆樣頻率Fs為1024Hz,釆樣點(diǎn)數(shù)為1024點(diǎn)，則可以分辨到1Hz。024Hz的釆樣率釆樣1024點(diǎn)，剛好是1秒，也就是說(shuō)，釆樣1秒時(shí)間的信號(hào)并做FFT,則結(jié)果可以分析到1Hz,如果釆樣2秒時(shí)間的信號(hào)并做FFT,則結(jié)果可以分析到0.5Hzo如果要提高頻率分辨力，則必須增加采樣點(diǎn)數(shù)，也即采樣時(shí)間。頻率分辨率和采樣時(shí)間是倒數(shù)關(guān)系。2?采樣出的點(diǎn)經(jīng)過(guò)FFT運(yùn)算后，各自點(diǎn)對(duì)應(yīng)的頻率是多少？答：每一個(gè)點(diǎn)就對(duì)應(yīng)著一個(gè)頻率點(diǎn)。這個(gè)點(diǎn)的模值，就是該頻率值下的幅度特性。具體跟原始信號(hào)的幅度有什么關(guān)系呢？假設(shè)原始信號(hào)的峰值為A,那么FFT的結(jié)果的每個(gè)點(diǎn)（除了第一個(gè)點(diǎn)直流分量之外）的模值就是A的N/2倍。而第一個(gè)點(diǎn)就是直流分量，它的模值就是直流分量的N倍。

3?信號(hào)的相位怎么通過(guò)FFT計(jì)算出的點(diǎn)推算出來(lái)？答：假設(shè)FFT之后某點(diǎn)n用復(fù)數(shù)a+bi表示，那么這個(gè)復(fù)數(shù)的模就是An二根號(hào)a*a+b*b,相位就是Pn=atan2(b,a)。(atan2)0幵七求相位為什么差90度menggxingg|分類；工程技術(shù)科學(xué)|瀏覽429次 2012-06-060分李到；S分李到；S舉報(bào)舉報(bào)▼ 2012-06-061L丄提問(wèn)者采納正弦和余弦之間的相位正好差90S，可能是這個(gè)原因；如果你的信號(hào)模型是正弦，你所用的方法求得的是余弦型的，那么就正好差90度了，沒(méi)問(wèn)題的,據(jù)我說(shuō)知，F(xiàn)FT求得的相位應(yīng)該是余弦型的atan2返回給定的X及Y坐標(biāo)值的反正切值。反正切的角度值等于X軸正方向與通過(guò)原點(diǎn)和給定坐標(biāo)點(diǎn)(￥坐標(biāo),邂標(biāo))的射線之間笊夾角。結(jié)果以弧度表示芥介于叩i到pi之間(不包括叩i)?？偨Y(jié)：根據(jù)以上的結(jié)果，就可以計(jì)算出n點(diǎn)(nHl,且nUN/2)對(duì)應(yīng)的信號(hào)的表達(dá)式為：An/(N/2)*cos(2*pi*Fn*t+Pn),即2*An/N*cos(2*pi*Fn*t+Pn)。對(duì)于nh點(diǎn)的信號(hào)，是直流分量，幅度即為A1/N。由于FFT結(jié)果的對(duì)稱性，通常我們只使用前半部分的結(jié)果，即小于采樣頻率一半的結(jié)果。借用高人用matlab做的一個(gè)來(lái)理解：假設(shè)我們有一個(gè)信號(hào)，它含有2V的直流分量，頻率為50Hz、相位為-30度、幅度為3V的交流信號(hào)，頻率為75Hz、相位為90度、幅度為1.5V的交流信號(hào)。用數(shù)學(xué)表達(dá)式就是如下：S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)式中cos參數(shù)為弧度，所以-30度和90度要分別換算成弧度。我們以256Hz的釆樣率對(duì)這個(gè)信號(hào)進(jìn)行釆樣，總共釆樣256點(diǎn)。按照我們上面的分析，F(xiàn)n=(n-l)*Fs/N,我們可以知道，每?jī)蓚€(gè)點(diǎn)之間的間距就是1Hz,第n個(gè)點(diǎn)的頻率就是n-l。我們的信號(hào)有3個(gè)頻率：0Hz、50Hz、75Hz,應(yīng)該分別在第1個(gè)點(diǎn)、第51個(gè)點(diǎn)、第76個(gè)點(diǎn)上出現(xiàn)峰值，其它各點(diǎn)應(yīng)該接近0。實(shí)際情況如何呢？很明顯，1點(diǎn)、51點(diǎn)、76點(diǎn)的值都比較大，它附近的點(diǎn)值都很小，可以認(rèn)為是0,即在那些頻率點(diǎn)上的信號(hào)幅度為0。接著，我們來(lái)計(jì)算各點(diǎn)的幅度值。分別計(jì)算這三個(gè)點(diǎn)的模值，結(jié)果如下：1點(diǎn)：51251點(diǎn)：38476點(diǎn)：192按照公式，可以計(jì)算出直流分量為：512/N二512/256二2；50Hz信號(hào)的幅度為:384/(N/2)=384/(256/2)=3：75Hz信號(hào)的幅度為192/(N/2)=192/(256/2)=1.5O可見(jiàn),從頻譜分析出來(lái)的幅度是正確的。然后再來(lái)計(jì)算相位信息。直流信號(hào)沒(méi)有相位可言，不用管它。先計(jì)算50Hz信號(hào)的相位,atan2(-192,332.55)=-0.5236,結(jié)果是弧度，換算為角度就是180*(-0.5236)/pi二-30.0001。再計(jì)算75Hz信號(hào)的相位,atan2(192,3.4315E-12)=1.5708弧度,換算成角度就是180*1.5708/pi二90.0002?？偨Y(jié)：假設(shè)釆樣頻率為Fs,采樣點(diǎn)數(shù)為N,做FFT之后，某一點(diǎn)n(n從1開(kāi)始)表示的頻率為：Fn=(n-l)*Fs/N；該點(diǎn)的模值除以N/2就是對(duì)應(yīng)該頻率下的信號(hào)的幅度(對(duì)于直流信號(hào)是除以N)；該點(diǎn)的相位即是對(duì)應(yīng)該頻率下的信號(hào)的相位。相位的計(jì)算可用函數(shù)atan2(b,a)計(jì)算。atan2(b,a)是求坐標(biāo)為(a,b)點(diǎn)的角度值，范圍從-pi到pi。要精確到xHz,則需要采樣長(zhǎng)度為1/x秒的信號(hào)，并做FFT。要提高頻率分辨率，就需要增加釆樣點(diǎn)數(shù)，這在一些實(shí)際的應(yīng)用中是不現(xiàn)實(shí)的，需要在較短的時(shí)間內(nèi)完成分析。解決這個(gè)問(wèn)題的方法有頻率細(xì)分法，比較簡(jiǎn)單的方法是采樣比較短時(shí)間的信號(hào)，然后在后面補(bǔ)充一定數(shù)量的0,使其長(zhǎng)度達(dá)到需要的點(diǎn)數(shù)，再做FFT,這在一定程度上能夠提高頻率分辨力。我寫(xiě)的另外一種C語(yǔ)言的版本用VC6.0運(yùn)行結(jié)果：(具體程序給我留言)模擬