【正文】 k(n) || mpym *AR7+,T3,AC0 。AC0+=b2k*wk(n2) macm *(AR3+T1),*AR7+,AC0 。AC0+=b0k*wk(n1) macm *AR3+,*AR7+,AC0 。AC0+=b1k*wk(n) filter_loop mov rnd(hi(AC02)),*AR1+ 。按 Q15格式存放結(jié)果sample_loop popboth XAR7 ?；謴?fù) AR7 pop T3 ?；謴?fù) T3 popm ST3_55 。恢復(fù) ST1,ST2,ST3 popm ST2_55 popm ST1_55 ret .end */圖 711 IIR濾波器系數(shù)和信號(hào)緩沖區(qū)配置 快速傅里葉變換 FFTn FFT算法原理n 庫(kù)利一圖基算法n FFT算法的實(shí)現(xiàn) FFT算法原理l 快速傅里葉變換（ FFT）是離散傅里葉變換（ DFT）的一種快速算法。通過(guò) FFT算法， DFT的計(jì)算量大大減少，運(yùn)算時(shí)間縮短 1~2個(gè)數(shù)量級(jí)。l DFT的變換公式為 其中 為旋轉(zhuǎn)因子。l FFT之所以減少運(yùn)算量，主要是利用了旋轉(zhuǎn)因子的以下 3點(diǎn)特性：u對(duì)稱性u(píng)周期性u(píng)可約性利用這些特性可以使 DFT運(yùn)算中有些項(xiàng)進(jìn)行合并，將長(zhǎng)序列的DFT分解為短序列的 DFT。DFT從算法上分為按時(shí)間抽選（ DIT）和按頻率抽選 (DIF）。基 2的 DIT又被稱為庫(kù)利一圖基算法。基 2的 DIF又稱為桑德 —圖基算法。 庫(kù)利一圖基算法n 信號(hào)流圖n 比特反轉(zhuǎn)n 蝶形運(yùn)算1. 信號(hào)流圖2．比特反轉(zhuǎn)l 圖 711的輸入信號(hào)的順序是按照比特反轉(zhuǎn)排列的，輸出序列是按照自然順序的。比特反轉(zhuǎn)就是將序列下標(biāo)用二進(jìn)制表示，然后將二進(jìn)制數(shù)按照相反的方向排列，即得到這個(gè)序列的實(shí)際位置。l 按照自然排序的時(shí)域信號(hào)數(shù)據(jù)是 x(0)、 x(1)、 x(2)、x(3)、 x(4)、 x(5)、 x(6)、 x(7)，其序號(hào)寫(xiě)成二進(jìn)制數(shù)分別為 000b、 001b、 010b、 011b、 100b、101b、 110b、 111b，將這些二進(jìn)制數(shù)前后倒轉(zhuǎn)，即得到進(jìn)行 FFT前數(shù)據(jù)所對(duì)應(yīng)的實(shí)際二進(jìn)制數(shù)地址：000b、 100b、 010b、 110b、 001b、 101b、 011b、 111b，對(duì)應(yīng)的十進(jìn)制數(shù)是： 0、 7。序號(hào)為 3的存儲(chǔ)單元，按照自然排序應(yīng)該存放x(3)，但由于 FFT計(jì)算規(guī)則的要求，現(xiàn)在應(yīng)該存放 x(6)。 l 基 2DIT FFT算法，共由 M級(jí)構(gòu)成，每級(jí)計(jì)算由 N/2個(gè)蝶形運(yùn)算構(gòu)成。l 基本運(yùn)算單元為以下所謂蝶形運(yùn)算：l 蝶形運(yùn)算中上下兩個(gè)節(jié)點(diǎn) p、 q的間距為： FFT算法的實(shí)現(xiàn)l 為了敘述簡(jiǎn)單，本書(shū)給出采用 C語(yǔ)言編寫(xiě)的 FFT程序，相應(yīng)的匯編程序請(qǐng)讀者自行完成。（ 1）主程序： include include /*包含浮點(diǎn)復(fù)數(shù)結(jié)構(gòu)體定義頭文件*/extern void bit_rev(plex *,unsigned int)。/*位反轉(zhuǎn)函數(shù)聲明 */extern void fft(plex *,unsigned int,plex *,unsigned int)。 /* fft函數(shù)聲明 */例 713，基 2DIT FFT算法的 C語(yǔ)言實(shí)現(xiàn)。extern void generator(int *,unsigned int)define N 128 /* FFT的數(shù)據(jù)個(gè)數(shù) */define M 7 /* M=log2(N) */define PI plex X[N]。 /* 說(shuō)明輸入信號(hào)數(shù)組 ,為復(fù)數(shù) */plex W[M]。 /* 說(shuō)明旋轉(zhuǎn)因子數(shù)組 e^(j2PI/N),為復(fù)數(shù) */ plex temp。 /* 說(shuō)明臨時(shí)復(fù)數(shù)變量 */float xin[N]。float spectrum[N]。 /* 說(shuō)明功率譜信號(hào)數(shù)組 ,為實(shí)數(shù) */float re1[N],im1[N]。 /* 說(shuō)明臨時(shí)變量數(shù)組 ,為實(shí)數(shù) */void main(){ unsigned int i,j,L,LE,LE1。 /* *//* 產(chǎn)生旋轉(zhuǎn)因子表 */ for (L=1。 L=M。 L++) { LE=1L。 /* 子 FFT中的點(diǎn)數(shù) LE=2^L */ LE1=LE1。 /* 子 FFT中的蝶形運(yùn)算數(shù)目 */ W[L1].re = cos(PI/LE1)。 W[L1].im = sin(PI/LE1)。 } /* */ generator(xin,N)。 for (。) { for (i=0。 iN。 i++) { /* 構(gòu)造輸入信號(hào)樣本 */ X[i].re =xin[i]。 X[i].im = 0。 /* 復(fù)制到參考緩沖器 */ re1[i] = X[i].re。 im1[i] = X[i].im。 } /* 啟動(dòng) FFT */ bit_rev(X,M)。 /* 以倒位次序排列 X[] */ fft(X,M,W,1)。 /* 執(zhí)行 FFT */ /* 計(jì)算功率譜，驗(yàn)證 FFT結(jié)果 */ for (i=0。 iN。 i++) { = X[i].re*X[i].re。 = X[i].im*X[i].im。 spectrum[i] = ( + )*4。 } }}（ 2）浮點(diǎn)復(fù)數(shù)基 2 DIT FFT函數(shù)： include void fft(plex *X,unsigned int M,plex *W,unsigned int SCALE){ plex temp。 /* 復(fù)變量臨時(shí)存儲(chǔ)器 */ plex U。 /* 旋轉(zhuǎn)因子 W^k */ unsigned int i,j。 unsigned int id。 /* 蝶形運(yùn)算中下位節(jié)點(diǎn)的序號(hào) */ unsigned int N=1M。 /* FFT 的點(diǎn)數(shù) */ unsigned int L。 /* FFT 的級(jí)序號(hào) */ unsigned int LE。 /* L級(jí)子 FFT的點(diǎn)數(shù) */ unsigned int LE1。 /* L級(jí)子 FFT蝶形運(yùn)算的個(gè)數(shù) */ float scale。 scale = 。 for(L=1。L=M。L++) { LE=1L。 LE1=LE1。 = 。 = 0.。 for(j=0。jLE1。j++) { for(i=j。iN。i+=LE) /*進(jìn)行蝶形計(jì)算 */ { id=i+LE1。 =(X[id].re*[id].im*)*scale。 =(X[id].im*+X[id].re*)*scale。 X[id].re=X[i].re*。 X[id].im=X[i].im*。 X[i].re=X[i].re*scale+。 X[i].im=X[i].im*scale+。 } /*遞推計(jì)算 W^k*/ =*W[L1].*W[L1].im。 =*W[L1].im+*W[L1].re。 =。 } }} （ 3）位反轉(zhuǎn)函數(shù)： include void bit_rev(plex *X,unsigned int M){ plex temp。 unsigned int i,j,k。 unsigned int N=1M。 /* FFT 的點(diǎn)數(shù) */ unsigned int N2=N1。 for (j=0,i=1。 iN1。 i++) { k=N2。 while(k=j){j=k。 k=1。 } j+=k。 if(ij) { temp=X[j]。 X[j]=X[i]。X[i]=temp。}}}（ 4）信號(hào)發(fā)生器函數(shù)： /* 該程序用于產(chǎn)生一組信號(hào)樣本 */include define PI= define Fs=8000 。采樣頻率設(shè)為 8000Hzdefine T=1/Fs 。采樣時(shí)間為 define f1=500 。信號(hào)源頻率 1取為 500Hzdefine a1= 。信號(hào)源幅度 1取為 define w1=2*PI*f1*T void generator(int *x,unsigned int N){ unsigned int i。 for(i=0。iN。i++) { X[i]=a1*cos((float)w1*i)。}}（ 5）復(fù)數(shù)結(jié)構(gòu)定義頭文件： struct cmpx{float re。float im。}。typedef struct cmpx plex。