数字信号处理 吴镇扬-CHP3快速傅立叶变换

第3章

傅立叶变换及其快速算法

之

快速傅里叶变换概述

快速傅里叶变换(FFT)是求解离散傅里叶变换(DFT)的快速算法。问题的提出

直接计算N点DFT需要的计算量是多少？

计算一个X(k)需要N次复数乘法和N一1次复数加法。算出全部N点X(k)共需N2次复数乘法和N(N一1)次复数加法.

总运算量近似地正比于N2

。当N值很大（如2-D图像处理），运算量将非常庞大；同时，对于实时性很强的信号处理来说，必将对计算速度有十分苛刻的要求。为此，需要改进对DFT的计算方法，以减少总的运算次数。概述在正交矩阵中，虽然有N2个元素，但只有N个不同的值，且有些取值特别简单，主要由于旋转因子具有如下的特点：对称性周期性下面以四点DFT为例来说明快速算法的思路。如何充分利用这些关系？概述概述交换矩阵第二列和第三列得从上面的结果可以看出,利用对称性和周期性，求四点DFT只需要一次复数乘法，称为Coolkey-Tukey算法。概述算法分类：N为2的整次幂：按基数分为基-2FFT算法、基-4FFT算法、混合基FFT算法、分裂基FFT算法;当N不是2的整次幂:典型的有Winograd算法.按抽取方法分：时间抽取(Decimation－in－Time，简称DIT)；频率抽取(Decimation－in－Frequency，简称DIF)

概述3.3.1时间抽取（DIT）基2FFT算法

为了将大点数的DFT分解为小点数的DFT运算，要求序列的长度N为N＝2M(M为正整数)。该情况下的变换称为基2FFT。核心思想是N点DFTN/2点

DFTN/4点

DFT2点

DFT

1个2个4个N/2个问题是如何分最有效？可以对时间变量分（DIT)，也可对频率变量分(DIF）3.3.1时间抽取（DIT）基2FFT算法基本思路：从时域将N点序列x(n)按奇偶项分解为两组，分别计算两组N/2点DFT，然后再合成一个N点DFT，按此方法继续下去，直到2点DFT，从而减少运算量。算法具体步骤：一、算法的推导1、序列x(n)按奇偶项分解为两组，将DFT运算也相应分为两组则3.3.1时间抽取（DIT）基2FFT算法2、两个N/2点的DFT合成一个N点DFT问题：A(k)，B(k)都只有N/2个点，怎样得到X(k)的后N/2点？利用周期性和对称性得4.2时间抽取（DIT）基2FFT算法3.3.13.3.1时间抽取（DIT）基2FFT算法4.2时间抽取（DIT）基2FFT算法3、继续分解（一直分解到两点DFT变换）A(K)和B(K)仍是高复合数(N／2)的DFT，我们可按上述方法继续以分解。令r＝2l，r＝2l十1，l＝0，1，…，N／4-1，则A(K)和B(K)可分别表示为3.3.1时间抽取（DIT）基2FFT算法4.2时间抽取（DIT）基2FFT算法3.3.1时间抽取（DIT）基2FFT算法令则4.2时间抽取（DIT）基2FFT算法3.3.1时间抽取（DIT）基2FFT算法同理，令则按此方法一直分解下去直到2点DFT，当N=8时，如下：4.2时间抽取（DIT）基2FFT算法3.3.1时间抽取（DIT）基2FFT算法3.3.1时间抽取（DIT）基2FFT算法下面通过讨论寻找FFT的一般规律。二、算法的讨论1、“级”的概念在分解过程中，每分一次，称为一级运算。因为M=log2N，所以N点DFT可以分解为M级，按抽取算法的信号流图中来定义，从左到右分别称为0级、1级到M-1级。3.3.1时间抽取（DIT）基2FFT算法2、蝶形单元在算法的信号流图中，第m级存在这种运算，这种结构几何形状像蝴蝶，称为蝶形单元p、q是参于本蝶形单元运算的上、下节点的序号。由于第m级序号的两点只参于这一个蝶形单元的运算，其输出在第m十l级。且这一蝶形单元也不再涉及别的点。由于这一特点，在计算机编程时，我们可将蝶形单元的输出仍放在输入数组中，这一特点称为“同址运算”。4.2时间抽取（DIT）基2FFT算法3.3.1时间抽取（DIT）基2FFT算法3.3.1时间抽取（DIT）基2FFT算法

由于一级都含有N/2个蝶形单元，每个蝶形单元需要1次复数乘法和两次复数加法，因此完成log2N级共需要的复数乘法和加法分别为

直接计算DFT时所需的复乘数与复加数都是与N2成正比的。所以采用FFT算法使运算量大大减少。显然，N值愈大，节省的运算量愈多。3.3.1时间抽取（DIT）基2FFT算法3、“组”的概念在分解过程中，每一级的N/2个蝶形单元可以分成若干组，每一组具有相同的结构和W因子分布。第m级可分成N/2m+1组。例：N=8=23，分3级。第一级的分组及Wr因子如下：m=0级,分成四组：因子为m=1级,分成二组,因子为m=2级,分成一组,因子为3.3.1时间抽取（DIT）基2FFT算法4、Wr因子的分布由上分析可知结论：每由后向前（m由M-1-->0级）推进一级，则此系数为后级系数中偶数序号的那一半。3.3.1时间抽取（DIT）基2FFT算法5、码位倒置在FFT算法中，输出的频谱依照正常次序排列，但输入的序列x(n)是按奇偶分开的，分开的规律，以N=8为例，按如下方法进行排序（1）、将x(n)的序号写成二进制

x(000),x(001),…,x(110)

，x(111)。（2）将二进制的码进行翻转，得

x(000),x(100),…,x(011)

，

x(111)。（3）将二进制的翻转码转换为对应的十进制

x(0),x(4),…,x(3)，x(7)。这就是按奇偶抽取得到的顺序。3.3.1时间抽取（DIT）基2FFT算法说明：①在上述的基2FFT算法中，由于每一步分解都是按输入序列x(n)在时域上的次序是属于偶数还是奇数来抽取的，所以称为“按时间抽取法”或“时间抽取”。

②上述的基2FFT算法中，抽取也可在频域进行，引出频率抽取（DIF）基2FFT算法。3.3.2频率抽取（DIF）基2FFT算法

设输入序列长度为N=2M(M为正整数)，频率抽取法将输入序列不是按奇、偶分组，而是按前后对半分开，这样可将N点DFT写成前后两部分;将该序列的频域的输出序列X(k)(也是N点序列，按其频域顺序的奇偶分解为越来越短的子序列，称为基2按频率抽取的FFT算法。也称为Sander-Tukey算法。3.3.2频率抽取（DIF）基2FFT算法算法分析

现将输入x(n)按n的顺序分前后两部分:前半子序列x(n),0≤n≤N/2-1;后半子序列x(n+N/2),0≤n≤N/2-1;例：N=8时，前半序列为：x(0),x(1),x(2),x(3);后半序列为：x(4),x(5),x(6),x(7);考虑N点的DFT,由DFT定义得3.3.2频率抽取（DIF）基2FFT算法3.3.2频率抽取（DIF）基2FFT算法按k的奇偶将X(k)分成奇偶两部分,k=2r和k=2r+1,考虑k为偶数情况令3.3.2频率抽取（DIF）基2FFT算法考虑k为奇数情况令3.3.2频率抽取（DIF）基2FFT算法结论

一个N点的DFT被分解为两个N/2点;与时间抽取法的推演过程一样，由于N=2M,因此,N/2仍为偶数，所以可以将N/2点DFT的输出X(k)再分为偶数组和奇数组，这样就将一个N/2点的DFT分成两个N/4点DFT的输入，也是将N/2点的DFT的输入上、下对半分后通过蝶形运算而形成，直至最后为2点DFT。8点DIF基2FFT算法流图3.3.2频率抽取（DIF）基2FFT算法3.3.2频率抽取（DIF）基2FFT算法3.3.2频率抽取（DIF）基2FFT算法DIT与DIF的相同之处：（1）DIF与DIT两种算法均为原位运算。（2）DIF与DIT运算量相同。DIT与DIF的不同之处：（1）DIF与DIT两种算法结构倒过来。DIF为输入顺序，输出乱序。运算完毕再运行“二进制倒读”程序。DIT为输入乱序，输出顺序。先运行“二进制倒读”程序，再进行求DFT。（2）DIF与DIT根本区别：在于蝶形结不同。DIT的复数相乘出现在减法之前。DIF的复数相乘出现在减法之后。3.3.3N为组合数的FFT

以上讨论的都是以2为基数的FFT算法，即N=2M，这种情况实际上使用得最多。

优点：程序简单，效率高，使用方便。实际应用时，有限长序列的长度N很大程度上由人为因素确定，因此多数场合可取N=2M，从而直接使用以2为基数的FFT算法。如N不能人为确定,N的数值也不是以2为基数的整数次方,处理方法有两种:①补零:将x(n)补零,使N=2M.

例如N=30,补上x(30)=x(31)=0两点,使N=32=25,这样可直接采用以2为基数M=5的FFT程序。有限长度序列补零后并不影响其频谱X(ej)，只是频谱的采样点数增加了,上例中由30点增加到32点，所以在许多场合这种处理是可接受的。3.3.3N为组合数的FFT②如要求准确的N点DFT值,可采用任意数为基数的FFT算法,

其计算效率低于以2为基数FFT算法。如N为复合数，可分解为两个整数p与q的乘积,像前面以2为基数时一样,FFT的基本思想是将DFT的运算尽量分小,因此,在N=pq情况下,也希望将N点的DFT分解为p个q点DFT或q个p点DFT,以减少计算量。步骤：分别为0,1,…，Q-1;分别为0,1,…，P-1。3.3.3N为组合数的FFT

N点DFT可以重新写成为考虑到

再令令(1)先将x(n)通过x(n1Q+n0)改写成x(n1,n0)。因为

Q=4,n1=0,1,2,n0=0,1,2,3,故输入是按自然顺序的，即x(0,0)=x(0)x(0,1)=x(1)x(0,2)=x(2)x(0,3)=x(3)x(1,0)=x(4)x(1,1)=x(5)x(1,2)=x(6)x(1,3)=x(7)x(2,0)=x(8)x(2,1)=x(9)x(2,2)=x(10)x(2,3)=x(11)以P=3，Q=4,N=12为例

(2)求Ｑ个Ｐ点的DFT（３）X1(k0,n0)乘以得到X1′(k0,n0)。（４）求P个Ｑ点的DFT，参变量是k0（５）将X2(k0,k1)通过X(k0+k1P)恢复为X(k)3.3.3N为组合数的FFTN=12为组合数时的FFT3.3.3N为组合数的FFT（１）求Ｑ个P点DFT需要ＱP2次复数乘法和Ｑ·P·(P-1)次复数加法；（２）乘Ｎ个W因子需要Ｎ次复数乘法；（３）求P个Ｑ点DFT需要PQ2次复数乘法和P·Ｑ（Ｑ-1)次复数加法。

总的复数乘法量:QP2+N+PQ2=N(P+Q+1)；总的复数加法量:Q·P(P-1)+P·Q·(Q-1)=N(P+Q-2)例:N=23*29=667,N2=444889,N(P+Q+1)=35351当组合数

N=P1P2P3…Pm中所有的Pi均为4时，就是基四FFT算法

以上提出FFT算法，可以很快地求出全部DFT值。即求出有限长序列x(n)的z变换X(z)在单位园上N个等间隔抽样点zk处的抽样值。它要求N为高度复合数。即N可以分解成一些因子的乘积。例N=2L

实际上：（1）也许对其它围线上z变换取样发生兴趣。（2）只需要计算单位圆上某一段的频谱,即M不等于N。如窄带信号，希望在窄带频率内频率抽样能够非常密集，提高分辨率，带外则不考虑。（3）若N是大素数时，不能加以分解，又如何有效计算这种序列DFT。例N=311，若用基2则须补N=28=512点，要补211个零点。3.3.4线性调频Z变换3.3.4线性调频Z变换问题提出为了提高DFT的灵活性，须用新的方法。线性调频z变换(CZT)就是适用这种更为一般情况下，由x(n)求X(zk)的快速变换。CZT

来自于雷达专业的专用词汇。Z变换采用螺线抽样,可计算单位圆上任一段曲线的Z变换，适用于更一般情况下（M不等于N）由x(n)求X(zr)的快速算法,达到频域细化的目的，这种变换称为线性调频Z变换(简称CZT)。

为适应z可以沿平面内更一般的路径取值，我们沿z平面上的一段螺线作等分角的抽样，则z的取样点Zr可表示为：

已知N点序列x(n)，0≤n≤N-1,其z变换为其中M：表示欲分析的复频谱的点数。M不一定等于N。A，W都为任意复数,令

3.3.4线性调频Z变换一、CZT的定义3.3.4线性调频Z变换上式即为CZT的定义.现在讨论A0,W0,θ0,φ0的含义:为输出M点的变换域值.r=0时的A0ejθ0是CZT的起点,随着r的变化,r0,r1,…,rM-1构成CZT的变化路径，对于M-1点其极坐标为3.3.4线性调频Z变换3.3.4线性调频Z变换CZT在现Z平面上的变换路径是一条螺旋线。（1）A为起始样点位置起点半径，大于1时，表示螺旋线在单位圆外，反之，在单位圆内。起点半相角。（2）当W0>1,螺旋线内旋，反之外旋。（3）当A0=W0=1时，CZT的变换路径为单位圆上的一段弧，起点为P，终点为Q，且M不一定等于N。3.3.4线性调频Z变换3.3.4线性调频Z变换

我们希望做的是频谱分析，因此考虑A0=W0=1时，在单位圆上CZT，且M不一定等于N。3.3.4线性调频Z变换3.3.4线性调频Z变换CZT的线性滤波计算步骤3.3.4线性调频Z变换二、CZT的计算方法分析：从上面的推导过程可以看出，计算CZT关键是计算一个线性卷积其中，g(n)应为N点序列，h(n)应为偶对称的无限长序列，考虑到输出M点序列，h(n)的实际长度应为M点。因此，可用DFT来实现两者的卷积，具体步骤如下：3.3.4线性调频Z变换（1）计算并设置g(n)3.3.4线性调频Z变换（2）计算并设置h(n)将h(n)设置成长度为L的序列，考虑到其为偶对称序列，且取M点（输出M点序列），取如下图所示3.3.4线性调频Z变换3.3.4线性调频Z变换（3）计算h‘(n)和g’(n)的DFT，得到L点序列H‘(k)和G’(k)。（4）令Y‘(k)=H‘(k)G’(k)（乘积）后作Y‘(k)

的IDFT（反变换）得到时域输出序列y(r)

。（5）取y(r)的前M点，并乘以W-r2/2,则得最后的输出X(zr)，即

与标准FFT算法相比，CZT算法有以下特点：（1）输入序列长度N及输出序列长度M不需要相等，且N及M不必是高度合成数，二者均可为素数。（2）Zk的角间隔是任意的，说明其频率分辨率也是任意可控的，角间隔小，分辨率高，反之，分辨率低。（3）周线不必是z平面上的圆，在语音分析中螺旋周线具有某些优点。（4）由于起始点z0可任意选定，因此可以从任意频率上开始对输入数据进行窄带高分辨率的分析。总之，CZT算法具有很大的灵活性。3.3.4线性调频Z变换2.4FFT应用中的几个问题1)IDFT的运算方法IDFT：DFT：以上所讨论的FFT算法可用于IDFT运算——简称为IFFT比较IDFT的定义式：IDFT与DFT的差别：

1）把DFT中的每一个系数改为，

2）再乘以常数1/N，

第二种方法，完全不需要改动FFT程序，而是直接利用它作IFFT。考虑到故

IFFT计算分三步：①将X（k）取共轭（虚部乘以-1）

②对直接作FFT③对FFT的结果取共轭并乘以1/N，得x（n）。2）实数序列的FFT

以上讨论的FFT算法都是复数运算,包括序列x(n)也认为是复数,但大多数场合,信号是实数序列,任何实数都可看成虚部为零的复数,例如,求某实信号x(n)的复谱,可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用FFT求其离散傅里叶变换。这种作法很不经济,因为把实序列变成复序列,存储器要增加一倍,且计算机运行时,即使虚部为零,也要进行涉及虚部的运算,浪费了运算量。合理的解决方法是利用复数据FFT对实数据进行有效计算,下面介绍两种方法。

（1）用

一个N点FFT同时计算两个N点实序列的DFT

设x

(n)、y

(n)是彼此独立的两个N点实序列,且

X

(k)=DFT[x

(n)]，Y

(k)=DFT[y(n)]

则X

(k)、Y(k)可通过一次FFT运算同时获得。首先将x

(n)、y(n)分别当作一复序列的实部及虚部,令

g(n)=x

(n)+jy(n)通过FFT运算可获得g(n)的DFT值利用离散傅里叶变换的共轭对称性通过g(n)的FFT运算结果G(k),由上式也可得到Y(k)的值。

（2）用一个N点的FFT运算获得一个2N点实序列的DFT

设x(n)是2N点的实序列,现人为地将x(n)分为偶数组x1

(n)和奇数组x2

(n)

x1

(n)=x(2n)n=0,1,…,N-1

x2

(n)=x(2n+1)

n=0,1,…,N-1然后将x1

(n)及x2

(n)组成一个复序列：

y(n)=x1

(n)+jx2

(n)

通过N点FFT运算可得到：Y(k)=X1(k)+jX2(k)，N点根据前面的讨论，得到

为求2N点x(n)所对应X(k)，需求出X(k)与X1

(k)、X2

(k)的关系：

而1）由x1(n)及x2(n)组成复序列，经FFT运算求得Y(k)，2）利用共轭对称性求出X1(k)、X2(k)，3）最后利用上式求出X（k）,达到用一个N点的FFT计算一个2N点的实序列的DFT的目的。

X(k)=X1(k)+W2NkX2(k)所以

3)线性卷积的FFT算法

线性卷积是求离散系统响应的主要方法之一,许多重要应用都建立在这一理论基础上,如卷积滤波等。以前曾讨论了用圆周卷积计算线性卷积的方法归纳如下:

将长为N2的序列x(n)延长到L,补L-N2个零，将长为N1的序列h(n)延长到L,补L-N1个零，如果L≥N1+N2-1,则圆周卷积与线性卷积相等,此时,可用FFT计算线性卷积，方法如下:

a.计算X(k)=FFT[x(n)]b.求H(k)=FFT[h(n)]c.求Y(k)=H(k)X(k)k=0--L-1d.求y(n)=IFFT[Y(k)]n=0--L-1

可见，只要进行二次FFT,一次IFFT就可完成线性卷积计算。计算表明,L>32时,上述计算线性卷积的方法比直接计算线卷积有明显的优越性,因此,也称上述循环卷积方法为快速卷积法。

上述结论适用于x(n)、h(n)两序列长度比较接近或相等的情况,如果x(n)、h(n)长度相差较多,例如,h(n)为某滤波器的单位脉冲响应,长度有限,用来处理一个很长的输入信号x(n)，或者处理一个连续不断的信号，按上述方法,有三个问题：（1）h(n)要补许多零再进行计算,计算量有很大的浪费，或者根本不能实现。（2）系统的存储量要求极高。（3）带来了很大的系统延迟。为了克服上述三个问题，保持快速卷积法的优越性,可将x(n)分为许多段,每段的长度与h(n)接近,处理方法有两种：(1)

重叠相加法——由分段卷积的各段相加构成总的卷积输出h(n)x(n)则输入序列可表为：于是输出可分解为：

其中假定xi(n)表示x（n）序列的第i段：

1）先对h（n）及xi（n）补零，补到具有N点长度，N=N1+N2-1。一般选N=2M。由于yi（n）的长度为N1，而xi（n）的长度为N2，因此相邻两段yi（n）序列必然有N-N2=N1-1点发生重叠。2）用基2FFT计算yi（n）=xi（n）*h（n）。3）重叠部分相加构成最后的输出序列。计算步骤：a.

事先准备好滤波器参数H（k）=DFT[h(n)]，N点b.

用N点FFT计算Xi（k）=DFT[xi(n)]c.

Yi(k)=Xi(k)H(k)d.

用N点IFFT求yi（n）=IDFT[Yi(k)]e.

将重叠部分相加(2)重叠保留

这种方法和第一种方法稍有不同，即将上面分段序列中补零的部分不是补零，而是保留原来的输入序列值，这时，如利用DFT实现h(n)和xi(n)的循环卷积，则每段卷积结果中有N1-1个点不等于线性卷积值需舍去。

重叠保留法与重叠相加法的计算量差不多，但省去了重叠相加法最后的相加运算。

y0(n)中的[N1-1,L-1]点对应于线性卷积

x(n)*h(n)中的[0,N2-1]点y1(n)中的[N1-1,L-1]点对应于线性卷积

x(n)*h(n)中的[N2,2N2-1]点y2(n)中的[N1-1,L-1]点对应于线性卷积

x(n)*h(n)中的[2N2,3N2-1]点依此类推

，并将yi(n)拼接起来构成y

(n)

4）用FFT计算相关函数

相关的概念很重要，互相关运算广泛应用于信号分析与统计分析，如通过相关函数峰值的检测测量两个信号的时延差等。两个长为N的实离散时间序列x（n）与y（n）的互相关函数定义为

则可以证明，rxy(τ)的离散傅里叶变换为

Rxy（k）=X*(k)Y(k)

其中X(k)=DFT[x(n)],Y(k)=DFT[y(n)],Rxy(k)=DFT[rxy(τ)]，0≤k≤N-1互相关函数定义为

x(n)及y(n)的卷积公式相比较，我们可以得到相关和卷积的时域关系：

因

得

证毕。

当x（n）=y（n）时，得到x（n）的自相关函数为：

维纳——辛钦定理：自相关函数与信号功率谱互为傅立叶变换对。

上面的推导表明，互相关和自相关函数的计算可利用FFT实现。由于离散傅里叶变换隐含着周期性，所以用FFT计算离散相关函数也是对周期序列而言的。直接做N点FFT相当于对两个N点序列x(n)、y(n)作周期延拓，作相关后再取主值（类似圆周卷积）。实际一般要求的是两个有限长序列的线性相关，为避免混淆，需采用与循环卷积求线性卷积相类似的方法，先将序列延长补0后再用上述方法。利用FFT求两个有限长序列的线性相关：(5)

对R（k）作IFFT，取后N-1项，得取前N项，得（1）

设x(n)和y(n)的长均为N，求线性相关；（2）为了使两个有限长序列的线性相关可用其循环相关代替而不产生混淆，选择周期L≥2N-1，且L=2m，以使用FFT，将x(n)，y(n)补零至长为L。（3）

用FFT计算X(k)，Y(k)（k=0，1…，L-1）（4）

R（k）=X*（k）y（k）

x=[13-112331];

y=[21-1120-13];

k=length(x);

xk=fft(x,2*k);

yk=fft(y,2*k);

rm=real(ifft(conj(xk).*yk));

rm=[rm(k+2:2*k)rm(1:k)];

m=(-k+1):(k-1);

stem(m,rm)

xlabel('m');ylabel('幅度');

-8-6-4-202468-6-4-2024681012m幅度两个序列的自相关函数例8例9-20-1001020-50050100150200250300m-20-1001020-50050100150200250300m幅度幅度

(a)(b)延迟序列的互相关函数(a)和自相关函数(b)

5）用FFT计算二维离散的傅里叶变换

二维信号有图象信号、时空信号、时频信号等。二维离散傅里叶变换可用于处理二维离散信号。二维离散傅里叶变换的定义为：

二维离散傅里叶变换可通过两次一维离散傅里叶变换来实现:1）作一维N点DFT（对每个m做一次，共M次）

k=0,1,…,N-1，m=0,1,…,M-12）作M点的DFT（对每个k做一次，共N次）

k=0,1,…,N-1，