深圳pcb抄板实数FFT算法及其C语言的设计实现案例

    鉴于目前在许多嵌入式系统中要用到FFT运算，如以DSP为核心的交流采样系统、频谱分析、相关分析等。本人结合自己的实际开发经验，研究了实数的FFT算法并给出具体的C语言函数，读者可以直接应用于自己的系统中。

    首先分析实数FFT算法的推导过程，然后给出一种具体实现FFT算法的C语言程序，电路板克隆可以直接应用于需要FFT运算的单片机或DSP等嵌入式系统中。

    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代替。