基于DSP的FFT算法實現(xiàn)

艾 紅，常青青，鄧大偉

（北京信息科技大學(xué) 自動化學(xué)院，北京 100192）

0 引言

快速傅立葉變換（FFT）在雷達(dá)、通信、電子對抗和電力系統(tǒng)等領(lǐng)域有廣泛應(yīng)用，特別是在電力系統(tǒng)的諧波檢測中，F(xiàn)FT幾乎是唯一可行的檢測方法。通常提高FFT運(yùn)算速度有兩種途徑：改進(jìn)FFT算法本身和改進(jìn)運(yùn)算工具。現(xiàn)階段提高FFT算法本身非常困難，一般方法致力于改進(jìn)運(yùn)算工具。數(shù)字信號處理器DSP 是一種可編程的高性能處理器。文中充分利用TMS320F2812 DSP強(qiáng)大的數(shù)據(jù)處理能力，實現(xiàn)了FFT運(yùn)算，并提高了運(yùn)算速度。

1 系統(tǒng)硬件結(jié)構(gòu)

系統(tǒng)設(shè)計以TMS320F2812處理器為核心，輔以外圍電路構(gòu)成。DSP負(fù)責(zé)模擬輸入信號數(shù)據(jù)采集以及FFT算法實現(xiàn)。外圍電路包括電源轉(zhuǎn)換電路，時鐘電路，復(fù)位電路以及外部RAM等。系統(tǒng)的硬件整體結(jié)構(gòu)如圖1所示[1]。對信號進(jìn)行FFT變換，首先要對模擬信號采樣將其轉(zhuǎn)換為數(shù)字信號。輸入的連續(xù)模擬信號經(jīng)信號調(diào)理電路后輸出到DSP的ADC模擬輸入通道，經(jīng)過ADC數(shù)據(jù)采集，模數(shù)轉(zhuǎn)換的結(jié)果存放于ADC結(jié)果寄存器中。信號調(diào)理電路主要是為了信號的抗混疊濾波以及電路阻抗的匹配 。信號調(diào)理電路對輸入信號進(jìn)行調(diào)理處理，包括信號的濾波、跟隨輸出以及信號的穩(wěn)定。DSP對采集的數(shù)字信號進(jìn)行FFT運(yùn)算處理，同時對運(yùn)算結(jié)果進(jìn)行相應(yīng)的數(shù)據(jù)顯示和數(shù)據(jù)存儲。

圖1 系統(tǒng)的硬件整體結(jié)構(gòu)圖

2 FFT算法原理

FFT是離散傅立葉變換（DFT）的快速運(yùn)算，是數(shù)字信號處理的基礎(chǔ)。因為有些信號在時域很難看出特性，使用FFT將其變換到頻域，就會很容易看出其特性。DFT算法的基本公式為[3]：

FFT算法是不斷地把長序列的DFT分解成幾個短序列的DFT，并利用的周期性和對稱性減少DFT的運(yùn)算次數(shù)。設(shè)序列x(n)的長度為N(N=2M,M為任意正整數(shù))，按n的奇偶把x(n)分解成兩個N/2點(diǎn)的子序列：

由此可見，若將任何一偶數(shù)點(diǎn)序列按下標(biāo)的奇偶性分成兩個子序列，則原序列的DFT可由兩子序列的DFT線性組合得到。運(yùn)算流圖如圖2所示。

圖2 運(yùn)算流圖

圖3 蝶式運(yùn)算流圖

其中，A和B的距離稱為翅尖距。這種方法和直接進(jìn)行DFT計算相比較，運(yùn)算量減少一半。按照這種分解運(yùn)算的思想，將X1(k)和X2(k)繼續(xù)向下分解，直到最后變?yōu)橐稽c(diǎn)序列，此時的運(yùn)算量大大減小。以8點(diǎn)序列X(n)的FFT運(yùn)算為例：

蝶式運(yùn)算流圖如圖3所示。比較FFT和DFT的運(yùn)算量。假設(shè)序列點(diǎn)數(shù)為N，且有N=2L,當(dāng)使用FFT進(jìn)行運(yùn)算時，計算過程中只有蝶式運(yùn)算，它的復(fù)數(shù)乘法運(yùn)算量為，復(fù)數(shù)加法運(yùn)算量為；直接進(jìn)行DFT運(yùn)算時，復(fù)數(shù)乘法和復(fù)數(shù)加法的運(yùn)算量均為N2。由此可見，F(xiàn)FT運(yùn)算確實大大減小了DFT的運(yùn)算量。

由圖3可以看出，若想得到順序正確的頻域序列X(k)，必須對時域序列進(jìn)行重新排序。這里先說明比特逆序的概念：設(shè)存儲地址m,其二進(jìn)制數(shù)（設(shè)m=2L）為m=(m1m2m3),若有=（m3m2m1），則稱為m的L位比特逆序列。若將圖4中的輸入時域序列下標(biāo)轉(zhuǎn)換成二進(jìn)制，依次是：000,100,010,110,001,101,011,111；而原時域序列下標(biāo)轉(zhuǎn)換成二進(jìn)制依次是:000,001,010,011,100,101,110,111，將兩者對比可發(fā)現(xiàn)，輸入時域序列下標(biāo)就是原時域序列下標(biāo)的比特逆序。所以FFT蝶式運(yùn)算的第一步就是對時域序列進(jìn)行比特排序。

3 系統(tǒng)軟件實現(xiàn)以及結(jié)果分析

3.1 FFT程序?qū)崿F(xiàn)

FFT算法主要包括比特逆序、蝶式權(quán)值的計算、各級的蝶式運(yùn)算以及FFT序列的輸出。程序中涉及到復(fù)數(shù)的運(yùn)算，由于計算機(jī)無法處理復(fù)數(shù)，故而將復(fù)數(shù)拆為實部和虛部，只需計算出實部和虛部的數(shù)值即可。程序設(shè)計將蝶式權(quán)w[i]拆分為實部pr[i]和虛部pi[i]，將參與各級蝶式運(yùn)算的序列拆分為實部fr和虛部fi。復(fù)數(shù)進(jìn)行加法運(yùn)算時需將兩復(fù)數(shù)的實部與實部、虛部與虛部對應(yīng)相加即可，例如w[i]+w[i+1],所得的和的實部為pr[i]+pr[i+1]，虛部為pi[i]+pi[i+1]；而復(fù)數(shù)的乘法，相對來說有些復(fù)雜，不再是簡單的虛部和實部各自相乘,例如w[i]×w,所得的乘積的實部為pr[i]×pr+pi[i]×pi，虛部為:(pr[i]+pi[i])×(pr+pi)-( pr[i]×pr+pi[i]×pi)。該程序執(zhí)行的是128個點(diǎn)的FFT運(yùn)算，共有7級蝶式運(yùn)算，每個蝶式權(quán)都可看作是的n次方，所以在一個循環(huán)內(nèi)即可計算出所有的蝶式權(quán)。每級蝶式運(yùn)算計算出新的序列值代替該級的輸入序列值，這樣避免了開辟新的存儲空間，有效地節(jié)省了存儲空間。

FFT算法程序設(shè)計如下所示：

Uint16 ConversionCount;

Uint16 px[128];

Uint16 pz[128];

void kfft(pr,pi,n,k,fr, fi ,l,il)

{ Uint16 n,k,l,il;

double pr[],pi[],fr[], fi [];

int it,m,is,i,j,nv,l0;

double p,q,s,vr,vi,poddr,poddi;

for (it=0; it＜=n-1; it++)

{ m=it;

is=0;

for (i=0; i＜=k-1; i++)

{ j=m/2;

is=2*is+(m-2*j);

m=j;}

fr[it]=pr[is];

fi[it]=pi[is]; //此循環(huán)為比特逆序的實現(xiàn)，pr、pi是原序列的實部與虛部

} //fr、 fi 是排完序后的序列

pr[0]=1.0; pi[0]=0.0;

p=6.283185306/(1.0*n);

pr[1]=cos(p);

if (l!=0) pi[1]=-pi[1];

for (i=2; i＜=n-1; i++)

{ p=pr[i-1]*pr[1];

q=pi[i-1]*pi[1];

s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);

pr[i]=p-q;

pi[i]=s-p-q;} //此循環(huán)為權(quán)值運(yùn)算實現(xiàn)

for (it=0; it＜=n-2; it=it+2) //it是序列下標(biāo)，

在第一級蝶式運(yùn)算中

//相鄰的兩個運(yùn)算蝶下標(biāo)相差2

{ vr=fr[it]; vi= fi [it];

fr[it]=vr+fr[it+1]; //第一級蝶式運(yùn)算中相鄰兩個數(shù)值進(jìn)行加減運(yùn)算

fi [it]=vi+ fi [it+1];

fr[it+1]=vr-fr[it+1];

fi [it+1]=vi- fi [it+1];} //此循環(huán)為第一級蝶式運(yùn)算，該級中的蝶式權(quán)均為1

m=n/2; nv=2;

for (l0=k-2; l0＞=0; l0--) //除去第一級，還有六級蝶式運(yùn)算

{ m=m/2; nv=2*nv; // 同一級相鄰的兩個蝶式運(yùn)算中數(shù)值的下標(biāo)相差為nv

for (it=0; it＜=(m-1)*nv; it=it+nv) //it仍為數(shù)值的下標(biāo)

for(j=0;j＜=(nv/2)-1;j++)

{ p=pr[m*j]*fr[it+j+nv/2]; // nv/2為蝶式運(yùn)算的翅間距

q=pi[m*j]*f i [it+j+nv/2];

s=pr[m*j]+pi[m*j];

s=s*(fr[it+j+nv/2]+f i [it+j+nv/2]);

poddr=p-q; poddi=s-p-q;

fr[it+j+nv/2]=fr[it+j]-poddr;

fi [it+j+nv/2]=f i [it+j]-poddi;

r[it+j]=fr[it+j]+poddr;

fi [it+j]=f i [it+j]+poddi; //計算出的新序列值代替輸入序列值,節(jié)省了存儲空間

}}} // fr、f i 分別為計算出的FFT序列的實部和虛部

3.2 運(yùn)行結(jié)果

用FFT的計算公式對單一頻率的正余弦信號進(jìn)行變換，可以得出單一頻率的信號在頻域表現(xiàn)為在其正負(fù)頻率點(diǎn)上的兩個脈沖。根據(jù)單頻信號的這一特性可驗證上述程序是否正確。將單一頻率的正弦波信號接入硬件系統(tǒng)，打開圖形觀察窗口，將觀察點(diǎn)數(shù)設(shè)置為128點(diǎn)，即可觀察到如圖4所示的波形。

圖4 時域正弦波的FFT運(yùn)行結(jié)果

圖5 改善后的FFT運(yùn)行結(jié)果

圖4中上半部分波形為輸入的時域正弦波，下半部分波形為經(jīng)過變換后的頻域波形。在圖4中，由于產(chǎn)生的模擬信號伴有噪音，導(dǎo)致信號頻率不單一，故而在信號的頻域中，除了信號對應(yīng)的頻率點(diǎn)有脈沖外，其它頻率點(diǎn)處也有較弱的脈沖。

對輸入信號的頻率進(jìn)行調(diào)整，使信號的組成頻率盡可能單一，可以得到更好的波形，如圖5所示。可以看到在信號對應(yīng)的頻率點(diǎn)有脈沖，其它頻率點(diǎn)的振幅為0，與理論推導(dǎo)結(jié)果一致，由此可驗證程序設(shè)計的正確性。

4 結(jié)論

FFT是聲學(xué)、圖像和信號處理等領(lǐng)域中一種重要的分析工具，文中闡述了硬件結(jié)構(gòu)圖和信號調(diào)理電路。詳細(xì)介紹了FFT算法原理，采用C語言編寫程序?qū)崿F(xiàn)了FFT算法。程序運(yùn)行結(jié)果表明TMS320F2812 DSP實現(xiàn)FFT運(yùn)算速度快，精度高[4]。

[1] 賈瑋,楊錄,張艷花.基于TMS320VC5416的FFT算法的實現(xiàn)[J].山西電子技術(shù),2009,2:11-13.

[2] Rright Hert. Rapid algorithms of digital signal processing[M].Beijing,Electronic Industrial Press,2002.

[3] 胡廣書.數(shù)字信號處理[M].北京:清華大學(xué)出版社,2003.

[4] 蘇奎峰,呂強(qiáng),耿慶鋒,陳圣儉.TMS320F2812原理與開發(fā)[M].北京:電子工業(yè)出版社,2005.