FFT的C语言算法实现

FFT计算的理论公式:将N的数值的FFT按照奇数和偶数项拆分为两个N/2的FFT,推出一个迭代关系

《FFT的C语言算法实现》

将离散DFT某点的傅里叶变换看做偶数位置和奇数位置离散点的傅里叶变换的合成,如上述的合成关系,W^k是二者的合成桥梁。

以八点时间抽取基2FFT为例,其butterfly结构图如下:

《FFT的C语言算法实现》

  在执行蝶形运算之前,先要将数据按码位倒置排序,这样在FFT的结果中才会按自然顺序排列。

《FFT的C语言算法实现》

把N点复输入序列或2N点的是序列两两组合为复数按照实部,虚部,实部,虚部……的顺序在存储器中存储,下面将它们的顺序按照码位倒置的顺序重新排列。

码位倒置的实现方法:(1)简单的利用按位与、或循环实现

                                                             (2)利用如下的公式推导出第I项的逆序数与第I-1项逆序数的关系,利用迭代的方法,提高程序的执行效率。

《FFT的C语言算法实现》

             例如(偶数) 01110000——>00001     110                             将该偶数减去1再求其逆序数

(奇数)01110000-1=01101111——>11110     110后面的三项相同,前k+1项不同,在这个奇数的逆序数基础上减去前K项的加权和在减去第K+1项的加权,即可以得到此偶数对应的逆序数

由此可以得到奇数项的逆序数等于前一项的逆序数加上2^(n-1),而偶数项的逆序数等于前一项的逆序数加上上述的修正项得到,实现程序如下:

《FFT的C语言算法实现》

码位倒置之后数据在存储器中的结构如下(16点实序列):

《FFT的C语言算法实现》

蝶形运算部分程序如下:

《FFT的C语言算法实现》

FFT 算法的基本思想分析

  我们知道N点FFT运算可以分成log2(N)级,每一级都有N/2个碟形,FFT的基本思想是用3层循环完成全部运算(N点FFT)。

  第一层循环:由于N=2^m需要m级计算,第一层循环对运算的级数进行控制。(stages)

  第二层循环:由于第L级有2^(L-1)个蝶形因子(乘数),第二层循环根据乘数进行控制,保证对于每一个蝶形因子第三层循环要执行一次,这样,第三层循环在第二层循环控制下,每一级要进行2^(L-1)次循环计算。(选择W)

  第三层循环:由于第L级共有N/2^L即2^(n-L)个羣,并且同一级内不同羣的乘数分布相同,当第二层循环确定某一乘数后,第三层循环要将本级中每个羣中具有这一乘数的蝶形计算一次,即第三层循环每执行完一次要进行N/2^L个碟形计算。(执行不同group中具有相同W的蝶形运算)

  可以得出结论:在每一级中,第三层循环完成N/2^L个碟形计算;第二层循环使第三层循环进行 2^(L-1)次,因此,第二层循环完成时,共进行2^(L-1) *N/2^L=N/2个碟形计算。实质是:第二、第三层循环完成了第L级的计算。

       其中的istep是具有相同W乘数旋转因子的子变换在存储器中的距离,这里实部、虚部相间排列,故此距离应增加一倍。W旋转因子采用迭代的方法实现(这里的程序是别人的,所以未作替换),当然也可以利用matlab产生后放到数组中,计算时查表获取,可参照TI的FFT库中的做法。

             最后将istep再付给mmax,由于istep时计算蝶形运算时两数据的间隔,当其大于点数N的时候,外层循环结束。个人觉得用stages更清晰些。

 

 

点赞