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