分类： C/C++
2011-07-28 19:56:31
快速傅立叶变换(FFT)
4．1引言
快速傅立叶变换（FFT）并不是一种新的变换，而是离散傅立叶变换（DFT）的一种快速算法。
DFT的计算在数字信号处理中非常有用。例如在FIR滤波器设计中会遇到从h(n)求H(k)或由H(k)计算h(n)，这就要计算DFT；信号的谱分析对通信、图像传输、雷达等都是很重要的，也要计算DFT。因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大。
自从1965年图基（J. W. Tukey）和库利（T. W. Coody）在《计算数学》（Math. Computation , Vol. 19, 1965）杂志上发表了著名的《机器计算傅立叶级数的一种算法》论文后，桑德(G. Sand)-图基等快速算法相继出现，又经人们进行改进，很快形成一套高效运算方法，这就是快速傅立叶变换简称FFT（Fast Fourier Transform）。这种算法使DFT的运算效率提高1~2个数量级。
4．2 基2 FFT算法
一、直接计算DFT的问题及改进的途径
设x(n)为N点有限长序列，其DFT正变换为
= ， k=0，1，…，N-1
其反变换(IDFT)
x(n)= ，n=0，1，…，N-1
二者的差别只在于 的指数符号不同，以及差一个常数乘因子1/N，因而下面我们只讨论DFT正变换的运算量，反变换的运算量是完全相同的。
考虑x(n)为复数序列的一般情况，每计算一个X(k)，需要N次复数乘法以及（N-1）次复数加法。因此，对所有N个k值，共需N2次复数乘法及
N（N-1）次复数加法运算。所以直接计算DFT，乘法次数和加法次数都是和N2成正比的，当N很大时，运算量是很可观的，因而需要改进对DFT的计算方法，以减少运算次数。
下面讨论减少运算工作量的途径。仔细观察DFT的运算就可看出，利用系数 以下固有特性，就可减小DFT的运算量：
（1） 的对称性 （ ）*=
（2） 的周期性 = =
（3） 的可约性 = =
由此可得： = = ， =-1， =- 。这样利用这些特性，可以将长序列的DFT分解为短序列的DFT。快速傅立叶变换算法正是基于这样的思路而发展起来的。它的算法基本上可以分成两大类：时域抽取法FFT（Decimation-In-Time FFT，简称DIT-FFT）和频域抽取法FFT（Decimation-In-Time FFT，简称DIF-FFT）。
二、时域抽取法基-2 FFT原理
先设序列点数为N=2M，M为整数。如果不满足这个条件，可以人为地加上若干零值点，使之达到这一要求。这种N为2的整数幂的FFT称基-2 FFT。
设输入序列长度为N=2M (M为正整数) ，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为按时间抽取(DIT )的FFT算法。也称Cooley - Tukey算法。
将N=2M的序列x(n)先按n的奇偶分成以下两组:
x(2r)=x1(r) r=0，1，…，N/2-1
x(2r+1)=x2(r)
则x(n)的DFT为：
＝DFT[x(n)]= ，k=0，1，…，N-1
= +
= +
= +
= +
= + ，k=0，1，…，N-1 （4.2.4）
式中 与 分别是x1(r)及x2(r)的N/2点DFT：
= = ，k=0，1，…，N/2-1
= = ，k=0，1，…，N/2-1
由（4.2.4）式可看出，一个N点DFT已分解成两个N/2点的DFT，它们按（4.2.4）式又组合成一个N点DFT。
现讨论 的周期性：
= = =
式中用了 =
同理可得： =
再考虑 的以下性质：
= =－
所以前半部分：
= + ，k=0，1，…，N/2-1
后半部分：
= +
= － ，（k=0，1，…，N/2-1）
这样，只要求出0到（N/2-1）区间的所有 和 值，即可求出0到（N-1）区间内的所有 值，这就大大节省了运算。
= + ，k=0，1，…， (4.2.7)
= － ，k=0，1，…， (4.2.8)
采用蝶形运算符号表示的图示法，可将上面讨论的分解过程表示于下图中。此图表示N=23=8的情况，其中输出值X（0）到X（3）是由（4.2.7）式给出的，而输出值X（4）到X（7）是由（4.2.8）给出。
既然如此，由于N=2M，因而N/2仍是偶数，可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。
x1(2r) = x3(r) r=0，1，…，N/4-1
x1(2r+1)=x4(r)
= +
= +
= + ，k=0，1，…，N/4-1
且 = - ，k=0，1，…，N/4-1
同理， 也可进行同样的分解：
= + ，k=0，1，…，N/4-1
且 = - ，k=0，1，…，N/4-1
其中， = =
= =
进一步具体化：N=8=23，
=
=x(0)+x(4) ，k=0，1
= x(0)+ x(4)
= x(0)+ x(4) = x(0)- x(4)
注意上式中 =
同理， = x(1)+ x(5)
= x(1)- x(5)
三、DIT-FFT算法与直接计算DFT运算量的比较
由上图得，当N=2M时，其运算流图有M级蝶形，每一级都有N/2个蝶形运算构成。每一个蝶形运算需要1次复数乘和2次复数加。所以每一级运算都需要N/2次复数乘和N次复数加。
M级运算总共需要的复数乘法次数为：
M级运算总共需要的复数加法次数为：
如直接计算DFT，复数乘法为 次，复数加法 次。
当N=1024时，复数乘之比： ，这样使运算效率提高200多倍。下图为FFT算法和直接DFT算法所需运算量与计算点数N的关系曲线：
显然，N越大，优越性就越明显。
四、按时间抽选的FFT算法的特点：
1、 原位运算
由图4.2.4可以看出，DIT-FFT的运算过程很有规律。N=2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据结点又同在一条水平线上，这就意味着计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元。这样，经过M级运算后，原来存放输入序列数据的N个存储单元中便依次存放 的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位计算。原位计算可节省大量内存，从而使设备成本降低。
2、 倒序规律
由图4.2.4看出，按原位计算时，FFT的输出 是按正常顺序排列在存储单元中，即按 ， ，…， 的顺序排列，但是这时输入x(n)却不是按自然顺序存储的，而是按x(0)，x(4)，…，x(7)的顺序存入存储单元，看起来好象是“混乱无序”的，实际上是有规律的，我们称之为倒序。
造成倒序的原因是输入x(n)按标号n的偶奇的不断分组而造成。由于N=2M，所以倒序数可用M位二进制数（nM-1nM-2…n0）2（当N=8=23时，二进制为三位）表示。第一次分组，标为n0。 n为偶数在上半部分，用n0=0表示，n为奇数在下半部分，用n0=1表示；第二次分组，标为n1。偶数部分再分为偶（0）奇（1），奇数部分再分为偶（0）奇（1）…。依次类推，直到M次分组，最后所得二进制倒序数如图示。
下表列出了N=8时以二进制数表示的顺序数和倒序数，由表显而易见，只要将顺序数（n2n1n0）的二进制位倒置，则得对应的二进制倒序值（n0n1n2）。
3、 倒序的实现
设原输入序列x(n)先按自然顺序存入数组A中。例如N=8，A(0)，A(1)，A(2)，A(3)，A(4)，A(5)，A(6)，A(7)中依次存放着x(0)，x(1)，x(2)，x(3)，x(4)，x(5)，x(6)，x(7)。
顺序数用I表示，I=1~N-2。倒序数用J表示，与I对应分别为4，2，6，1，5，3。当I=J时不需要交换，I
I=1时，对应的倒序数是4；I=2时，对应的倒序数是2…。倒序数从4到2到…的关系可从表4.2.1得到：每次最高位加1。（注意J用十进制数表示）。如果最高位为0，J直接加N/2，如果最高位为1，则要将最高位归0，次高位加1。但次高位加1时也要判断是否为1或0。程序框图如下图虚线框里所示：
4、 蝶形运算两个输入数据的“距离”
以图4.2.4的8点FFT为例，其输入是倒位序的，输出是自然顺序的。N=2M，共有第一级蝶形运算，第二级蝶形运算，…，第M级蝶形运算。用L表示第某级运算，每个蝶形的两个输入数据的“距离”为B=2L-1。
5、 旋转因子的变化规律
仍观察图4.2.4，每级都有N/2个蝶形，每个蝶形都要乘以因子 ，称其为旋转因子，p称为旋转因子的指数。第L级蝶形运算（从左向右数），共有2L-1个
L=1，第一级： = ，J=0
L=2，第二级： = ，J=0，1
L=3，第三级： = ，J=0，1，2，3
所以对于N=2M的一般情况，第L级的 为：
，J=0，1，…，2L-1-1
由于2L=2M2L-M=N2L-M ，所以 = = =
所以，第L级的 为： ，J=0，1，…，2L-1-1
例如：
L=1，第一级：2L-1=1，J=0，
L=2，第二级：2L-1=2，J=0，1， ，
L=3，第三级：2L-1=4，J=0，1，2，3， ， ， ，
6、蝶形运算规律
设序列x(n)经时域抽选（倒序）后，存入数组X中。如果蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶形运算可表示成如下形式：
式中 ；J=0，1，…，2L-1-1；
L=1，…，M。
DIT-FFT运算和程序框图如下：
同一旋转因子对应着间隔为2L点的2M-L个蝶形。
五、按频率抽选（DIF）的基2 FFT算法
设序列点数为N=2M ，M为整数。
＝ ，k=0，1，…，N-1
先把输入按n 的顺序分成前后两半：
＝ ，k=0，1，…，N-1
= +
= +
= ，k=0，1，…，N-1
= ，k=0，1，…，N-1
=(-1)k= 1，k为偶数
-1，k为奇数
将 分解成偶数组与奇数组，
当k取偶数即k=2r时（r=0，1，…，N/2-1），
= =
当k取奇数即k=2r+1时（r=0，1，…，N/2-1），
=
=
令 n=0，1，…，N/2-1
= ，r=0，1，…，N/2-1
= (4.2.16)
、 和 之间可用下图所示的蝶形运算符号表示：
（4.2.16）用下图表示，N=8。一次分解流图。
由于N=2M，N/2仍是偶数，继续将N/2点DFT分成偶数组和奇数组。图4.2.12表示N=8时二次分解运算流图。
最后完整的分解流图为下图：
这种算法是对 进行奇偶抽取分解的结果，所以称之为频域抽取法FFT。
DIF-FFT算法与DIT-FFT算法类似，不同的是DIF-FFT算法输入序列为自然顺序，而输出为倒序排列。另外，蝶形运算略有不同，DIT-FFT蝶形先乘后加，而DIF-FFT蝶形先加后乘。
上述两种FFT的算法流图形式不是唯一的。只要保证各节点所连支路及其传输系数不变，改变输入与输出点以及中间结点的排列顺序，就可以得到其他变形的FFT运算流图。
六、IDFT的高效算法
离散傅立叶反变换
x(n)= ，n=0，1，…，N-1
与离散傅立叶正变换（ = ， k=0，1，…，N-1）相比，只要将DFT中的系数 改变为 ，最后乘以1/N，就是IDFT的运算公式。流图输入为 ，输出为x(n)。因此原来的DIT-FFT改为IFFT后称为DIF-IFFT更合适；原来的DIF-FFT改为IFFT后称为DIT-IFFT更合适。
下图是由DIF-FFT运算流图改成的IFFT运算流图（DIT-IFFT）：
在实际中，有时为了防止运算过程中发生溢出，将1/N分配到每一级蝶形运算中。由于1/N= ，所以每级的每个蝶形输出支路均有一相乘因子1/2。如下图示：
如果希望直接调用FFT子程序计算IFFT，则可用下面的方法：
由于 x(n)=
∴ x*(n)=
x(n)=
= {DFT[ ]}*
这样可以先将 取共轭，然后直接调用FFT子程序进行DFT运算，最后取共轭并乘以1/N得到序列x(n)。
4．3 进一步减少运算量的措施
研究进一步减少运算量的途径，以程序的复杂度换取计算量的进一步提高
一、 多类蝶形单元运算
由图4.2.4已得出结论，N=2M点FFT共需要MN/2次复数乘法。
由 = ，J=0，1，…，2L-1-1得，当L=1时，只有一种旋转因子 ，所以第一级不需要乘法运算。当L=2时，共有两个旋转因子： =1和 =-j，因此第二级也不需要乘法运算。在DFT中，称其值为±1和±j的旋转因子为无关紧要的旋转因子。
综上所述，先除去第一、第二两级后，所需复数乘法次数
进一步考虑各级中的无关紧要旋转因子。当L=3时，有两个无关紧要的旋转因子 和 。因为同一旋转因子对应着2M-L=N/2L个蝶形运算，所以，第三级共有2N/23=N/4个蝶形不需要复数乘法运算。依次类推，当L≥3时，第L级的2个无关紧要的旋转因子减少复数乘法的次数为2N/2L=N/2L-1。这样，L=3至L=M共减少复数乘法次数为
因此
=
在基2 FFT程序中，若包含了所有旋转因子，则称该算法为一类蝶形单元运算；若去掉 =±1的旋转因子，则称之为二类蝶形单元运算；若再去掉 =±j的旋转因子，则称为三类蝶形单元运算；若再处理 = ，则称之为四类蝶形运算。我们将后三种运算称为多类蝶形单元运算。显然蝶形单元越多，编程就越复杂，但当N较大时，乘法运算的减少量是相当可观的。例如，N=4096时，三类蝶形单元运算的乘法次数为一类蝶形单元运算的75%。
二、 旋转因子的生成
在FFT运算中，旋转因子
= ，求余弦和正弦函数值的计算量很大，所以编程时，一种方法是在每级运算中直接产生，另一种方法是在FFT程序开始前预先计算好，存放在数组中，作为旋转因子表，在程序执行过程中直接查表得到。
三、 实序列的FFT算法
在实际工作中，数据x(n)一般都是实序列。如果直接按FFT运算流图计算，就是把x(n)看成一个虚部为零的复序列进行计算，这就增加了运算时间。处理这个问题有二种方法，一种是早期提出的用一个N点FFT计算N点实序列的FFT。第二种方法是用N/2点FFT计算一个N点实序列的DFT。
