数字信号处理第四章-快速傅里叶变换

1FFT:FastFourierTransform1965年，JamesW.Cooley和JohnW.Tukey在《计算数学》（《MathematicsofComputation》）上发表了“一种用机器计算复序列傅立叶级数的算法（AnalgorithmforthemachinecalculationofcomplexFourierseries）”论文。自此之后，新的算法不断涌现。一种是对N等于2的整数次幂的算法，如基2算法，基4算法。另一种是N不等于2的整数次幂的算法，例如Winagrad算法，素因子算法。4.1FFT：引言Dr.JamesW.CooleyWorkedatIBMWatsonResearchCenterinYorktownHeights,N.Y..AfterhisretirementfromIBMin1991,hejoinedtheElectricalEngineeringDepartmentattheUniversityofRhodeIsland.24.2FFT：直接计算DFT的问题及改进N点有限长序列x(n)的DFT变换对的定义为：

其中假设x(n)是复序列，同时X(k)一般也是复数。3复数乘法复数加法每一个X(k)NN–1N个X(k)(N点DFT)N2N(N–1)实数乘法实数加法一次复乘42一次复加2每一个X(k)4N2N+2(N–1)=2(2N–1)N个X(k)(N点DFT)4N22N(2N–1)4.2FFT：直接计算DFT的问题及改进复乘的加法次数复加的加法次数4如N=512、1024和8192时，DFT的乘法运算

N2=5122=218=262144（26万次）

N2=10242=220=1048576（105万次）

N2=81922=226=67108864（6千7百万次）对于大N，在实际中是不能接受的，无法“实时”应用DFT。一般地，在计算机中，一次加法比一次乘法所需时间要短;在DSP中，由于乘法用特殊的硬件电路专门完成，因此乘法和加法所需机器周期相同。Cooley与Turkey提出的FFT算法，大大减少了计算次数。如N=512时，FFT的乘法次数约为2000次，提高了约128倍，而且简化随N的增加而巨增，因而，用数值方法计算频谱得到实际应用。4.2FFT：直接计算DFT的问题及改进5FFT：WNkn

的性质x(0)x(2)x(1)x(3)-1-1-1-1W41X(0)X(1)X(2)X(3)6以4点DFT为例：直接计算需要：42=16

次复数乘。而按周期性及对称性，可以将DFT表示为：只需要

1

次复数乘FFT：WNkn

的性质7FFT算法分类:时间抽选法DIT:Decimation-In-Time频率抽选法DIF:Decimation-In-FrequencyFFT：基本思想基本思路：虽然存在不同的FFT方法，但其核心思想大致相同，即通过迭代，反复利用低点数的DFT完成高点数的DFT计算，以此达到降低运算量的目的。迭代：利用WNkn

的周期性、特殊点和对称性，合并DFT计算中很多重复的计算，达到降低运算量的目的。低点数：将傅氏变换DFT分解成相继小的DFT计算，即N变小，而计算量与N2

成正比。8算法原理设序列点数N=2M，M为整数。若不满足，则补零。

将序列x(n)按n的奇偶分成两组：N为2的整数幂的FFT算法称基-2FFT算法。即一组由偶数序号组成，另一组由奇数序号组成。4.3FFT：基2时间抽选法注意其长度为N/2910偶数取样点DFT为：

奇数取样点DFT为：

①k的整个范围为0～(N-1)，而X1(k)、X2(k)是由N/2个样点形成的DFT，x(2r)和x(2r＋1）的长度为N/2；②由这两个偶数和奇数N/2个时域样值可以计算出前N/2个DFT系数，也可以计算出后N/2个DFT系数。③问题：这前后N/2个DFT有无关系？k取N/2～(N-1)时，X1(k)、X2(k)、WN

情况如何？11观察则12因此：整个X(k)的计算，可以分解为前、后半部分的运算。而只要求出前一半，就可以由上式求出整个序列。同理而因此13上式表示为信号流程图：此信号流程图也称为蝶形流程图偶数点的DFT奇数点的DFT14以N=8为例，其信号流程图：偶数序列奇数序列15N/2仍为偶数，进一步分解：N/2N/416同理17以N=8为例，其信号流程图为N/2点DFTN/2点DFT18由于N=2M，这样逐级分解，直到2点DFT当N=8时，即分解到X3(k)，X4(k)，X5(k)，X6(k)，k=0,119每一次分解都是按时间域输入序列的奇偶性次序，分解为两个半长的子序列，所以称为按时间抽取法。仍以N=8为例：注：复数乘法5次，原来64次；复数加法为24次，原来56次。m=0级m=1级m=2级20

2m

点DFT的基2时间抽选计算过程可见：为什么引入FFT？出发点：考虑W因子的特点和性质，简化算法。DIT-FFT算法：时间域输入信号逐级分解为奇偶序列。上式中位序调整第0级蝶形复合第1级蝶形复合第m-1级蝶形复合x(n)……X(k)21例

使用基2时间抽选法FFT流图，计算下列数据的8点DFT：x=[4,-3,2,0,-1,-2,3,1]4-320-1-23142-13-30-214-123-3-201355-1-5-11-1355j-5-11j85+j-25-j-4-1+j-6-1-j85+j-25-j-4j√26jj√245+j+j√2-2+6j5-j+j√2125+j-j√2-2-6j5-j-j√222蝶形运算在FFT信号流图中每一级里节点都是成对存在的，计算时总是下节点的值乘以WNr，然后与上节点的值相加减，形成下一列两个节点值。这种计算的基本关系是：式中p,q是上下节点的序号。每一级中每对节点称为对偶节点，对偶节点的间距为：基2时间抽选法-蝶形运算23同址计算（同位特性）基2时间抽选法-同址计算1)不同级数的节点序号不变从蝶形运算可以看出第m列的xm(p),xm(q)经蝶形运算后，在第(m+1)列中的节点序号不变，即xm+1(p)中的p与xm+1(q)中的q仍分别是xm(p),xm(q)中的p和q值。242)蝶形运算是自成独立单元蝶形运算自成独立单元，即xm+1(p)、xm+1(q)只与xm(p)、xm(q)有关，而与其他节点的值无关，同时xm(p)、xm(q)也不参与另外的蝶形运算，后面的运算也不再用到前面的数据。因此，就可把计算结果xm+1(p),xm+1(q)放入计算前xm(p)、xm(q)的存储单元中，而不用再建新的存储单元—同址计算。同址计算所需存储量仅等于给定数据所需的存储量，这可大大节省存储单元。基2时间抽选法-同址计算25输入位序重排N点DFT分解为两个N/2点DFT输入序列按奇偶分组再分解再奇偶重排直到2点DFT。即FFT输出自然序列输入序列x(n)位序重排。基2时间抽选法-位序重排010011n0n1n226说明：首先把x(n)中的n表示为二进制。假定N=8，则 x(n)→x(n2n1n0)ni=0,1FFT的时间抽选法按n的奇偶分离而形成位置重排的原理如上图所示。由此可以看出，时间抽选法FFT的输入位序重排，输出分成上下两部分，输出结果为自然顺序。实际操作输入序列重排实际上就是完成二进制前后位序的相互交换：基2时间抽选法-位序重排27当N=2M时，共有M=log2N级蝶形；每级N/2个蝶形；每个蝶形有1次复数乘法，2次复数加法。复数乘法：复数加法：比较DFT/FFT

FFT：基2时间抽选法-运算量28FFT算法上面讨论的FFT算法假定N为2的整数次幂，即N=2M，是基2FFT算法。当N=Rv时，这样的算法叫做基RFFT算法，而当N=R1v1R2v2R3v3时，叫做混合基算法。当N是一个高度合数时，可得到最有效的算法。最受欢迎也最易编程的算法是基2FFT算法。对于不能分解的质数，或者当N不是2的整数次幂时，可以在信号的末尾补0，使其成为高度合数或2的整数次幂。MatlabFFT实现Matlab提供了内建的X=fft(x,N)

函数来计算矢量x的DFT。fft函数是机器码写成的，而不是以Matlab指令写成的，即不存在.m文件。因此，它的执行速度很快。采用混合算法。若N为2的幂，则得到高速的基2FFT算法；若N不是2的幂，则将N分解成质数，得到较慢的混合基FFT；若N为质数，则fft函数采用原始DFT算法。FFT：基2时间抽选法-Matlab实现29DFT定义表示为矩阵形式为：FFT：基2时间抽选法-矩阵表示30令：矩阵WN

为：DFT矩阵简单地表示为：特性FFT：基2时间抽选法-矩阵表示31FFT基2时间抽选法信号流图（N=8）-1W821×8点2×4点4×2点x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)-1-1-1-1-1-1-1-1-1-1-1W80W80W82W80W81W82W83C(0)C(1)D(0)D(1)E(0)E(1)F(0)F(1)A(0)A(1)A(2)A(3)B(0)B(1)B(2)B(3)m=0级m=1级m=2级位序重排FFT：基2时间抽选法-矩阵表示32FFT矩阵形式表示为：1)输入矢量x与P8

相乘，则只是输入的位序重排，没有乘法操作。FFT：基2时间抽选法-矩阵表示332)第0级运算：2点DFTFFT：基2时间抽选法-矩阵表示343)第1级运算：4点DFTFFT：基2时间抽选法-矩阵表示354)第2级运算：8点DFTFFT：基2时间抽选法-矩阵表示36FFT算法看成是DFT矩阵W8

的分解：这个因式分解对应于快速算法，因为矩阵F8(8)，F8(4)

，F8(2)和P8

的大部分元素为0，是稀疏矩阵。因此上述每个矩阵的乘法运算最多只需要8次复乘运算，而P8只是置换操作，没有乘法操作。不同的FFT算法对应于将DFT矩阵WN

分解成不同的稀疏矩阵。FFT：基2时间抽选法-矩阵表示37算法原理：根据时间-频率的对偶性时间抽选法：是把输入x(n)按奇偶分解成两个子序列，即N点x(n)序列→N/2点子序列，而输出X(k)是按自然顺序排列的。频率抽选法：是把输入x(n)按照前后对半分开，而不是奇偶数分开，而输出X(k)逐项分解成偶数点子序列和奇数点子序列。DFT变换表达式为：4.4FFT：基2频率抽选法(DIF-FFT)如果将输入x(n)按前后等分，即将求和分成两部分，范围分别为：38基2频率抽选法–算法原理39

按k的奇偶将X(k)分成两部分：40令:则X(2r)和X(2r+1)分别是x1(n)和x2(n)的N/2点DFT，记为X1(k)和X2(k)。41x1(0)x1(1)-1x1(2)x1(3)-1x2(0)x2(1)-1x2(2)x2(3)-1N/2点DFTN/2点DFTx(0)x(7)x(1)x(2)x(3)x(4)x(5)x(6)X1(0)=X(0)X2(0)=X(1)X1(1)=X(2)X1(2)=X(4)X1(3)=X(6)X2(1)=X(3)X2(2)=X(5)X2(3)=X(7)42FFT基2频率抽选法信号流图（N=8）逐级分解，直到2点DFT43基本蝶形运算不同DIT:先复乘后加减，W因子在上下节点都有体现DIF:先减后复乘，W因子仅体现在下节点FFT：基2时间和频率抽选法的异同44运算量相同同址计算FFT：基2时间和频率抽选法的异同DIF输出位序重排，DIT输入位序重排45DIT和DIF的基本蝶形互为信号流图转置DITDIFFFT：基2时间和频率抽选法的异同46输入倒位序，输出自然序（DIT）输入自然序，输出倒位序（DIF）输入输出均自然序各级具有相同几何形状输入倒位序，输出自然序输入自然序，输出倒位序FFT：基2抽选法-其它形式的流图47

FFT：基2抽选法-其它形式的流图48FFT：基2抽选法-其它形式的流图49FFT：基2抽选法-其它形式的流图50FFT：基2抽选法-其它形式的流图51FFT：基2抽选法-其它形式的流图52DFT和IDFT的定义：4.5FFT：IDFT快速算法(IFFT)DFT和IDFT的区别：①因子W的指数相差一个负号；②相差一个因子1/N。FFT算法中的分组方式、排序方式以及蝶形运算结构都可用于IFFT算法的设计，而这就是可依据现有的FFT算法直接得出IFFT算法的原因。53FFT和IFFT基本蝶形运算之间的关系设有序列x(n)，其DFT为X(k)，则IDFT为:

在FFT的时间抽选法中:

4.5FFT：IDFT快速算法(IFFT)对于IFFT算法，输入是X(k)和X(k+N/2)，输出是X1(k)和X2(k)。解上式可以得到X1(K)、X2(K)：X(k)X(k+N/2)X1(k)X2(k)-1548点DIT-IFFT算法-1x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)-1-1-1-1-1-1-1-1-10.5W8-00.5W8-20.5W8-00.5W8-10.5W8-20.5W8-3-1-10.5W8-00.5W8-20.50.50.50.50.50.50.50.50.50.50.50.5说明：1）分组过程是按时间序列x(n)的奇偶性在时域上展开的，故称此法为时间抽选算法DIT-IFFT；2）1/N的分解，N=2m

。4.5FFT：IDFT快速算法(IFFT)558点DIF-IFFT算法4.5FFT：IDFT快速算法(IFFT)56用FFT程序求IFFT的方法直接法：利用DIT、DIF的FFT程序，改变参量①把X(k)作为输入序列，而输出序列就是x(n)；②把因子WNkn

改为WN-kn

；③输入序列的每一个元素除以N。4.5FFT：IDFT快速算法(IFFT)57流程如下：

共轭法共轭FFT共轭乘1/N4.5FFT：IDFT快速算法(IFFT)58DFT不适用于：DFT是均匀分布在Z平面单位圆上N点处的频谱，如果取样点不均匀时，则很麻烦。只研究信号的某一频段，要求对该频段取样密集，提高分辨率(如窄带信号的频谱分析）；研究非单位圆上的取样值（如频谱峰值探测）；需要准确计算N点DFT，且N为大的素数；当x(n)是短时间序列时，则得到的频率分辨率2pi/N是很低的。提高频谱密度的办法：用补零的方法增加点数，但DFT的点数又大大增加，使计算工作量增大。CZT算法：对z变换采用螺旋线取样，chirp-z变换(线性调频z变换，Chirpz-transform)4.9线性调频z变换CZT及其快速算法59CZT的定义设x(n)为已知时间序列，其Z变换的形式为：式中复变量z为：这里s为拉普拉斯变量，是个实数CZT：定义60

按照上式计算Z[x(n)]必然是从z的实轴开始，为得到任意起始点和以螺旋线规律变化的z值，做如下假设：

其中，A、W为任意复数，θ0为A的起始角（第1个取样点（k=0）），A0为A起始半径，φ0

为在Z平面中相邻的zk

（即zk

和zk+1

）之间的夹角，W0

为任意正数值。M为要分析的点数，不一定等于N。61CZT表达式1）取样点沿螺旋线按角度间隔φ0分布，φ0>0时，取样轨迹逆时针旋转；φ0<0时，取样轨迹逆时针旋转。2）zk

周线是一条螺旋线：W0

表示螺旋线的伸展率；W0>1，随着k的增加，向内盘旋，朝向原点；W0<1，随着k的增加，向外盘旋；3)若W0=1，一段圆弧，若同时A0=1，则为单位圆一部分，此时CZT[x(n)]=DFT[x(n)]（M=N）。62CZT表示为线性卷积形式：利用布鲁斯坦(Bluestein)提出的等式

把在CZT[x(n)]中的Wnk

因子展开，得：CZT：快速算法计算量：直接计算M点的CZT，需要MN次复数乘法，M(N-1)次复数加法，需要快速算法。把上述运算转换为线性卷积形式，从而可以采用FFT算法，提高运算速度。63CZT：快速算法令则式中64CZT：快速算法g(n)和h(n)的取值范围：65CZT：快速算法h(n)x(n)g(n)X(zk)序列：可以想像为频率随时间（n）线性增长的复指数序列，在雷达中这类信号，称为线性调频信号（chirpSignal)，因此，这里z变换称为线性调频z变换。该信号的瞬时频率Ωi=2Ωt正比于时间t，因此，该信号的频率随时间线性增长。类比：66取值范围n：0～N-1，k：0～M-1因果序列g(n)的长度为N，而序列h(n)的长度为N+M-1，并位于区间内-(N-1)~(M-1)，所以卷积后的序列y(n)的长度为2N+M-2，且位于区间内–(N-1)～（N＋M－2）。实际上只要得到区间0～(M-1）内的线性卷积结果就能够完成CZT的计算。CZT：快速算法67用循环卷积代替线性卷积且不产生混迭失真的条件是：循环卷积的点数（周期）应大于或等于2N＋M－2，但CZT只需要前0～(M－1)值，对以后的其它值是否混迭并不关心。因此，可将循环卷积的点数缩减到最小为N+M-1。为了基2FFT运算，取L>=N+M-1，L=2m

。CZT：快速算法68g(n)、h(n)序列的构造为了进行循环卷积，把h(n)以周期L进行周期拓展，再取主值序列，从而得到圆周卷积的一个序列。从n=M开始补L－（N＋M－1）个零或任意值，补到n=L-N处。从n=L－N＋1到L－1，取h(n)的周期拓展序列h(L-n)。进行循环卷积的另一个序列g(n)只需要补零到L即可。69h(n)的构造：70CZT快速算法下面说明上述计算方案的具体实现，既然X(zk)可用线性卷积表示，我们就可以用DFT来计算线性卷积部分，只不过DFT换成FFT，从而得到CZT的快速算法。L点DFTH(k)h(n)构造长度M＋N-1长度Ng(n)补零L点DFTG(k)G(k)H(k)L=N+M-1L点IDFTy(k)7172线性卷积求解小结：时域直接求解

补N-N1个零x(n)N点DFT补N-N2个零h(n)N点DFTN点IDFTy(n)=x(n)*h(n)z变换法DFT法4.10FFT的应用：线性卷积求解73上面介绍的是两个有限长序列x1(n)、x2(n)的线性卷积。但有时其中某个序列会很长或无限长，若等长序列存储或输入完后再做卷积运算，将产生问题：

(1)存储量过大，运算量也太大；(2)等待输入的时间很长，引起不能忍受的延迟；(3)采用DFT/FFT快速算法的效率降低具体方法：重叠保留法、重叠相加法

解决此问题的思路：

把长序列分段，每一分段分别与短序列进行卷积——分段卷积。FFT的应用：分段卷积74误差分析在分析分段卷积之前，我们首先分析两个有限长度序列的循环卷积的误差情况。设x(n)为N1

点序列，h(n)为N2

点序列，由前面分析知道，线性卷积y(n)和循环卷积yC(n)关系为：DFT的应用：分段卷积一般地讲，循环卷积是线性卷积的一种混叠形式，即是线性卷积以周期N延拓后取主值周期。但当N≥N1+N2-1，没有混迭产生，此时线性卷积y(n)和循环卷积yN(n)相等。75为了用DFT计算线性卷积，必须选择适当的N，但实际中却可能做不到这一点，尤其是当一个序列的长度很大而存储空间有限时。当选择比所需要的小的N值进行循环卷积时，会引入误差。下面我们计算此误差的大小。首先N≥max(N1,N2)，假设此时，循环卷积yc(n)的取值范围为而线性卷积y(n)的非零范围为76有上式可知，误差为：77

因为max(N1,N2)≤N≤（N1＋N2－1），上面的求和式中只剩下对应于r=±1的项，其它的r值超出了max(N1,N2)≤N≤（N1＋N2－1）。因此，

一般来讲，如果x(n)和h(n)为因果序列，则y(n)也为因果序列，即：因此78它说明当

时，样本n处的误差与线性卷积中第n+N

个样本相等，即混叠部分的样值。（N1+N2–1）个样本后，线性卷积结果为零。这说明循环卷积中的头L个样本存在误差（混叠），混叠长度L等于线性卷积长度和循环卷积长度之差，而剩下的则为正确的线性卷积样本。

0循环卷积长度N(N1+N2-1)-Ny(n)线性卷积长度N1+N2-1循环卷积的错误样本（混叠部分）79此结论在用分段处理法实现长卷积时是非常有用的。错误样本数=线性卷积长度–

循环卷积长度N1N20假设N=N2(N1+N2-1)-N=N1-1e(n)=0e(n)=y(n+N)y(n)yC(n)h(n)x(n)N1+N2-1N1-1点有误差与线性卷积相同

假设L=N1-180例设x1(n)和x2(n)是两个四点序列：

x1(n)={1,2,2,1}，x2(n)={1,-1,-1,1}计算当N=6，5，4时的循环卷积，并在每种情况下，验证它们的误差。解：显然，线性卷积为7点x3(n)={1,1,-1,-2,-1,1,1}

当N=6，得到6点序列的循环卷积为：因此81当N=5时，得到5点序列的循环卷积为当N=4时，得到4点序列的循环卷积为82重叠保留法——对输入x(n)预处理基本思路假设把序列x(n)分成多段N点序列，一个系统的脉冲响应为M点序列，M<N。由上面的结论可知，输入分段序列和脉冲响应之间的N点循环卷积产生该段的输出序列，其中前M-1个样本不是正确的输出值。若将x(n)简单地分成互不重叠的各段，则所得的输出序列会有不正确样本区间存在。输入数据如何分段？结果如何处理？DFT的应用：重叠保留法MNy

(n)h(n)x(n)M-1点误差正确的点MNM-1点误差正确的点83为了解决这个问题，将x(n)分为长度N1

的分段，然后每段再往前多取（M-1）个样本，即第k段的最后M-1点保留下来作为第(k+1)段的前M-1点，各段长变为：

其中最前面的一段x0(n)，在其前面补上M-1个零，使其长度也为N。其中任一段xk(n)与h(n)进行N点卷积运算。

循环卷积：DFT的应用：重叠保留法相应的周期卷积 的周期为：N=N1+M-1线性卷积：84yk(n)的长度为：N1N’M-1M-1x0(n)nN重叠区域线性卷积循环卷积N1N’NM-1M-1n重叠区域线性卷积x1(n)错误样本数=线性卷积长度–

循环卷积长度

=[N1+2(M-1)]–[N1+(M-1)]=M-185将每一段所得到的循环卷积的前M-1个值去掉，保留后面的N1

个值，再将各段保留的N1个值前后拼接起来，就得到所要求的线性卷积y(n)。因此，循环卷积

yk’(n)的前M-1个值是线性卷积yk(n)前面的M-1个值与yk-1(n)后面M-1个值的混迭，是错误的样本。

循环卷积yk’(n)的后N1个值才与yk(n)中间的

N1个值相同。86

这种通过将与相重叠的部分舍去，即舍去循环卷积yk’(n)中前M-1项，保留后面的N1

个数值，来得到所求结果的方法称为重叠保留法，有的书中又称为重叠舍去法。8788例设x(n)=(n+1),0≤n≤9,h(n)={1,0,-1}，按N=6用重叠舍去法计算

y(n)=x(n)*h(n)解：由于M=3，必须使每一段与前一段重叠两个样本，x(n)为10点序列，需要在开头加(M-1)=2个零。因为N=