嵌入式系统中FFT算法研究

摘　要：首先分析实数fft算法的推导过程，然后给出一种具体实现fft算法的c语言程序，可以直接应用于需要fft运算的单片机或dsp等嵌入式系统中。
关键词：嵌入式系统 fft算法 单片机 dsp 目前国内有关数字信号处理的教材在讲解快速傅里叶变换（fft）时，都是以复数fft为重点，实数fft算法都是一笔带过，书中给出的具体实现程序多为basic或fortran程序并且多数不能真正运行。鉴于目前在许多嵌入式系统中要用到fft运算，如以dsp为核心的交流采样系统、频谱分析、相关分析等。本人结合自己的实际开发经验，研究了实数的fft算法并给出具体的c语言函数，读者可以直接应用于自己的系统中。 1 倒位序算法分析　　按时间抽取（dit）的fft算法通常将原始数据倒位序存储，最后按正常顺序输出结果x（0），x（1），...，x（k），...。假设一开始，数据在数组 float datar[128]中，我们将下标i表示为（b6b5b4b3b2b1b0）b，倒位序存放就是将原来第i个位置的元素存放到第（b0b1b2b3b4b5b6）b的位置上去.由于c语言的位操作能力很强，可以分别提取出b6、b5、b4、b3、b2、b1、b0，再重新组合成b0、b1、b2、b3、b4、b5、b6，即是倒位序的位置。程序段如下（假设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;　　大家可以对比教科书上的倒位序程序，会发现这种算法充分利用了c语言的位操作能力，非常容易理解而且位操作的速度很快。2 实数蝶形运算算法的推导　　我们首先看一下图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）。
设 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）];
继续分解得到下列两式:
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)的存储位置相同,所以经过（1）、（2）后，该位置上的值已经改变，而下面求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）]继续分解得到下列两式:
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代替。　　② 经过式（3）后， xr(k+b)的值已变化，而式（4）中要用到该位置上的上一级值，所以在执行式（3）前要先将上一级的值xr’(k+b)保存。　　③ 在编程时, xr(k)和 xr’(k), xi(k)和 xi’(k)使用同一个变量。
　　通过以上分析，我们只要将式（1）、（2）、（3）、（4）转换成c语言语句即可。要注意变量的中间保存，详见以下程序段。/* 蝶形运算程序段 ,datar[]存放实数部分,datai[]存放虚部*/
/* cos、sin函数做成表格，直接查表加快运算速度 */
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层循环完成全部运算(n点fft)。　　第一层循环：由于n=2m需要m级计算,第一层循环对运算的级数进行控制。　　第二层循环：由于第l级有2l-1个蝶形因子（乘数），第二层循环根据乘数进行控制，保证对于每一个蝶形因子第三层循环要执行一次，这样，第三层循环在第二层循环控制下，每一级要进行2l-1次循环计算。　　第三层循环：由于第l级共有