嵌入式系統(tǒng)中FFT算法研究

摘　要：首先分析實(shí)數(shù)fft算法的推導(dǎo)過(guò)程，然后給出一種具體實(shí)現(xiàn)fft算法的c語(yǔ)言程序，可以直接應(yīng)用于需要fft運(yùn)算的單片機(jī)或dsp等嵌入式系統(tǒng)中。
關(guān)鍵詞：嵌入式系統(tǒng) fft算法 單片機(jī) dsp 目前國(guó)內(nèi)有關(guān)數(shù)字信號(hào)處理的教材在講解快速傅里葉變換（fft）時(shí)，都是以復(fù)數(shù)fft為重點(diǎn)，實(shí)數(shù)fft算法都是一筆帶過(guò)，書(shū)中給出的具體實(shí)現(xiàn)程序多為basic或fortran程序并且多數(shù)不能真正運(yùn)行。鑒于目前在許多嵌入式系統(tǒng)中要用到fft運(yùn)算，如以dsp為核心的交流采樣系統(tǒng)、頻譜分析、相關(guān)分析等。本人結(jié)合自己的實(shí)際開(kāi)發(fā)經(jīng)驗(yàn)，研究了實(shí)數(shù)的fft算法并給出具體的c語(yǔ)言函數(shù)，讀者可以直接應(yīng)用于自己的系統(tǒng)中。 1 倒位序算法分析　　按時(shí)間抽?。╠it）的fft算法通常將原始數(shù)據(jù)倒位序存儲(chǔ)，最后按正常順序輸出結(jié)果x（0），x（1），...，x（k），...。假設(shè)一開(kāi)始，數(shù)據(jù)在數(shù)組 float datar[128]中，我們將下標(biāo)i表示為（b6b5b4b3b2b1b0）b，倒位序存放就是將原來(lái)第i個(gè)位置的元素存放到第（b0b1b2b3b4b5b6）b的位置上去.由于c語(yǔ)言的位操作能力很強(qiáng)，可以分別提取出b6、b5、b4、b3、b2、b1、b0，再重新組合成b0、b1、b2、b3、b4、b5、b6，即是倒位序的位置。程序段如下（假設(shè)128點(diǎn)fft）：
/* i為原始存放位置，最后得invert_pos為倒位序存放位置 */
　int b0=b1=b2=b3=b4=b5=6=0;
　b0=i&0x01; b1=(i/2)&0x01; b2=(i/4)&0x01;
　b3=(i/8)&0x01; b4=(i/16)&0x01; b5=(i/32)&0x01;
　b6=(i/64)&0x01; /*以上語(yǔ)句提取各比特的0、1值*/
　invert_pos=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;　　大家可以對(duì)比教科書(shū)上的倒位序程序，會(huì)發(fā)現(xiàn)這種算法充分利用了c語(yǔ)言的位操作能力，非常容易理解而且位操作的速度很快。2 實(shí)數(shù)蝶形運(yùn)算算法的推導(dǎo)　　我們首先看一下圖1所示的蝶形圖。蝶形公式：
x(k) = x’(k) + x’(k+b)w pn ,
x(k+b) = x’(k) - x’(k+b) w pn
其中w pn= cos（2πp/n）- jsin（2πp/n）。
設(shè) x(k+b) = xr(k+b) + jxi(k+b),
x(k) = xr(k) + jxi(k) ,
有:
xr(k)+jxi(k)= xr’(k)+jxi’(k)+[ xr’(k+b) + jxi’(k+b)]*[ cos（2πp/n）-jsin（2πp/n）];
繼續(xù)分解得到下列兩式:
xr(k)= xr’(k)+ xr’(k+b) cos（2πp/n）+ xi’(k+b) sin (2πp/n) （1）
xi(k)= xi’(k)-xr’(k+b) sin（2πp/n）+xi’(k+b)cos (2πp/n) （2） 　　需要注意的是: xr(k)、xr’(k)的存儲(chǔ)位置相同,所以經(jīng)過(guò)（1）、（2）后，該位置上的值已經(jīng)改變，而下面求x(k+b)要用到x’(k)，因此在編程時(shí)要注意保存xr’(k)和xi’(k)到tr和ti兩個(gè)臨時(shí)變量中?！　⊥? xr(k+b)+jxi(k+b)= xr’(k)+jxi’(k)- [ xr’(k+b)+jxi’(k+b)] *[ cos（2πp/n）-jsin（2πp/n）]繼續(xù)分解得到下列兩式:
xr(k+b)= xr’(k)-xr’(k+b) cos（2πp/n）- xi’(k+b) sin (2πp/n) （3）
xi(k+b)= xi’(k)+ xr’(k+b) sin（2πp/n）- xi’(k+b) cos (2πp/n) （4）
注意：
　　① 在編程時(shí), 式（3）、（4）中的xr’(k)和 xi’(k)分別用tr和ti代替?！　、?經(jīng)過(guò)式（3）后， xr(k+b)的值已變化，而式（4）中要用到該位置上的上一級(jí)值，所以在執(zhí)行式（3）前要先將上一級(jí)的值xr’(k+b)保存?！　、?在編程時(shí), xr(k)和 xr’(k), xi(k)和 xi’(k)使用同一個(gè)變量。
　　通過(guò)以上分析，我們只要將式（1）、（2）、（3）、（4）轉(zhuǎn)換成c語(yǔ)言語(yǔ)句即可。要注意變量的中間保存，詳見(jiàn)以下程序段。/* 蝶形運(yùn)算程序段 ,datar[]存放實(shí)數(shù)部分,datai[]存放虛部*/
/* cos、sin函數(shù)做成表格，直接查表加快運(yùn)算速度 */
tr=datar[k]; ti=datai[k]; temp=datar[k+b];/*保存變量，供后面語(yǔ)句使用*/
datar[k]=datar[k]+datar[k+b]*cos_tab[p]+datai[k+b]*sin_tab[p];
datai[k]=datai[k]-datar[k+b]*sin_tab[p]+datai[k+b]*cos_tab[p];
datar[k+b]=tr-datar[k+b]*cos_tab[p]-datai[k+b]*sin_tab[p];
datai[k+b]=ti+temp*sin_tab[p]-datai[k+b]*cos_tab[p];3 dit fft 算法的基本思想分析　　我們知道n點(diǎn)fft運(yùn)算可以分成logn2 級(jí)，每一級(jí)都有n/2個(gè)碟形。dit fft的基本思想是用3層循環(huán)完成全部運(yùn)算(n點(diǎn)fft)?！　〉谝粚友h(huán)：由于n=2m需要m級(jí)計(jì)算,第一層循環(huán)對(duì)運(yùn)算的級(jí)數(shù)進(jìn)行控制?！　〉诙友h(huán)：由于第l級(jí)有2l-1個(gè)蝶形因子（乘數(shù)），第二層循環(huán)根據(jù)乘數(shù)進(jìn)行控制，保證對(duì)于每一個(gè)蝶形因子第三層循環(huán)要執(zhí)行一次，這樣，第三層循環(huán)在第二層循環(huán)控制下，每一級(jí)要進(jìn)行2l-1次循環(huán)計(jì)算?！　〉谌龑友h(huán)：由于第l級(jí)共有