【转】数字信号处理教案 ---- 4、快速傅立叶变换(FFT)

dontium

【转】数字信号处理教案 ---- 4、快速傅立叶变换(FFT) [复制链接]

第四章 快速傅立叶变换(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

由（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

=(-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。

文章出处：http://wlsyzx.yzu.edu.cn/kcwz/szxhcl/kechenneirong/jiaoan/jiaoan4.htm