数字信号处理原理与实践 第4章

第四章快速傅里叶变换数字信号处理原理与实践(修订版)@NEU离散傅里叶变换（DFT）在数字信号处理中是最常用的运算之一，由于DFT计算比较繁琐，在相当长的时间内并没有得到真正的应用。直到1965年库利(J.W.Cooley)和图基(J.W.Tukey)提出了DFT的一种快速算法，又有桑德(G.Sande)和图基的快速算法相继出现，之后又出现了各种各样的DFT算法，这些方法统称为快速傅里叶变换（FFT）。4.1离散傅里叶变换存在的问题对于N点序列，其DFT为一般情况下，和都是复数，每计算一个值，需要作N次复数乘法和(N-1)次复数加法，求N点的值，需要N2次复数乘法，以及N(N-1)次复数加法。由于实现一次复数乘法需要四次实数乘法和两次实数加法，一次复数加法需要两次实数加法。

在实际运算过程中，运算量并没有如此之大，因为在DFT运算中包含着大量的重复运算，同时一些系数是不需要乘法运算的。4.2按时间抽取的基2FFT算法

利用的周期性，对称性，把长度为N点的DFT运算逐次分解为较短序列的DFT运算。这种方法叫时间抽取法(decimation-in-time,缩写为DIT)，简称时选法。4.2.1算法的推导

设N＝2M，M为正整数，N是2的M次方，如果不满足这个条件，可以人为的加上若干个零值点，达到这个要求。这种N为2的整数次幂的FFT称为基2FFT。将按照奇数和偶数分解为两个序列，将N点的DFT运算分解为两个N/2点的DFT运算，即令n=2r，n=2r+1，r=0，1，2，….N/2-1DFT运算分为对应的两组，于是有

k=0，1，…N/2-1k=0，1，…N/2-1k=0，1，…N/2-1令则后半部分的这样，N点的就可以用和来表示，分别为N/2个点的DFT。运算可以用蝶形信号流图表示。求N点的值，需要N2次复数乘法和N(N-1)次复数加法。对于N/2点DFT有(N/2)2=N2/4次复数乘法和(N/2)(N/2-1)次复数加法，因此，两个N/2点DFT有N2/2次复数乘法和N(N/2-1)次复数加法。把两个N/2的DFT合成一个N点DFT时，需要N/2个蝶型运算，包括N/2次复数乘法和N次复数加法。由于N>>1，总共的复数乘法为，复数加法次数为，通过这样的分解，运算量几乎可以减少到一半。可以按上述方法继续分解，将N/2点的DFT按奇偶部分分解为两个N/4点的DFT，进一步减少复数乘法的运算量。以N＝8为例，将一个8点的DFT分解为两个4点的DFT，其分解过程如图所示同样，再把每个4点的DFT分解为两个点的DFT，对于两点的DFT来说，已经不需要再分解。因此得到以下结果：当N＝8时，将一个N/2点的DFT分解成两个N/4个点的DFT的分解过程如图所示。这样一个8点的DFT就分解成为四个2点的DFT。如果N＝16，32，或者2的更高次幂，则可以按上述方法继续分解下去，直到只有两点DFT为止。4.2.2算法的讨论

用上面所述的方法，任何一个N＝2M点的DFT，都可以经过M次分解，最后分成两个点的DFT运算，如图所示，每分一次，称为一“级”运算，所以N点的DFT可以分为M级。由于每一级都含有N/2个蝶型单元，每一个蝶型单元又只需要一次复数乘法，两次复数加法，则完成一次时选FFT的总计算量为：复数乘法次数：Mf=(N/2)log2N=MN/2复数加法次数：Af=Nlog2N=MN当全部FFT完成后，变换后输出的序列依照正序排列，即存储单元中顺序放着X(0),X(1),…,X(N-1)。但输入却不是原来的自然顺序，而是x(0),x(4),x(2),…,x(N-1)是由于作奇偶分开所产生的。以N＝8为例，其自然序号是0，1，2，3，4，5，6，7。第一次按奇偶分开，得到两组N/2点的DFT，的序列号是0，2，4，6和1，3，5，7。对每一组再做奇偶分开，这时抽取后得到四组，每组的序号是0，4和2，6和1，5和3，7。这种从十进制看来很乱的输入顺序,实际上是按二进制“倒序位”排列的。仍然以N=8为例,x(0),x(1),…,x(7)对应是x(000),x(001),x(010),x(011),x(100),x(101),x(110),x(111),将二进制码翻转得到x(000),x(100),x(010),x(110),x(001),x(101),x(011),x(111),它们对应的十进制序号分别是x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)。4.3按频率抽取基2FFT算法

与按时间抽取相对应，常用的FFT算法还有基2频率抽取（decimation-in-frequency,缩写为DIF）FFT算法，它不是将按时序的奇偶分解，而是将频域的序号按奇偶分开，下面说明频率抽取的基2FFT算法。令序列点数N＝2M，M为整数，先将N点的DFT分成上下两个部分，即式中将X(k)分解为偶数组和奇数组,分别令k=2r,k=2r+1,r=0,1,2,…,N/2-1,则令得到两个N/2点DFT这样，一个N点的DFT被分解为两个N/2点的DFT。当N＝8时，上述过程如图所示：将一个N点DFT分解为4个N/4点DFT的过程如图所示。是一个完整的N＝8，按频率抽选的FFT结构如图所示。4.4运算量进一步减少的方法

虽然前面讨论的FFT算法，已经大大减少了运算量，但是还可以通过一些措施，进一步减少运算量。对于N＝2M，一共需要M级运算，每级N/2个蝶型单元，每个蝶型单元进行一次复数乘法，则总共需要MN/2次复数乘法，如果＝1，－1，j，-j，则与这些因子相乘时不作实际的复数乘法运算。4.5IDFT的快速计算方法IFFT

IDFT的快速计算方法是快速傅里叶反变换，用英文缩写IFFT来表示，无论时选法还是频选法的FFT算法均可用于IDFT运算。对于N点序列，其DFT运算公式为其IDFT运算公式为只要将DFT运算中的改为，在将结果乘以1/N，就可以用FFT的计算方法来实现IFFT。频域序列为输入，时域序列为输出，所以时选FFT对应时选的IFFT，频选FFT对应频选IFFT。对于系数1/N要做如下处理，把1/N分解为（1/2）M，在M的每一级运算都要乘以1/2，这样做的目的是为了避免在运算过程中出现溢出。其中N为IFFT变换的点数，N＝2M，且M为整数。IFFT的流图如图4.6分裂基FFT算法

分裂基FFT算法又称基2/4算法，或混合基算法，它和前面讲的基2算法有关，也和基4算法有关。基4算法比基2FFT算法更快，从理论上说，用较大的基数可以进一步减少运算次数，但程序过于复杂，所以一般不取基数大于8。1984年，法国的杜梅尔和霍尔曼首先提出了“分裂基”算法，把基2算法和基4算法结合起来，其基本思路是对偶数序号输出采用基2算法，对奇数数序号输出采用基4算法。由于分裂基算法是目前已知的对于N＝2M算法中具有最少的加法次数和乘法次数，并且具有非常好的结构，因此被认为是最好的FFT算法。后来的研究表明，该算法最接近理论上的所需乘法次数的最小值。首先介绍一下基4FFT算法。如果N＝16，通过上述推导可以把16点的DFT分成4个4点的DFT，其信号第一级流图如图所示。可以用同样的方法对X1(k)、X2(k)、X3(k)做第三级分解,把r为偶序号的x1(r)作N/4点的DFT,把r为奇序号的x1(r)作N8点的DFT。同样的,对X2(k)、X3(k)做同样的处理。以N=16为例,X(k)的第一级的分解可以得到4个分裂基,X1(k)的第二级分解可以得到2个分裂基,一个基-4的4点DFT和2个基-2的2点DFT,而X2(k)和X3(k)的第二级的分解分别是基-4的4点DFT,对N=16的分裂基FFT的示意图如图所示。4.7快速傅里叶变换的MATLAB程序实现

MATLAB为计算快速傅里叶变换，提供了一系列丰富的数学函数，主要有fft，Ifft，fft2，ifft2等。当所处理的数据的长度为2的幂时，采用基—2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2的整数次幂或者用添零方式来添补数据使其长度成为2的整数次幂。1．fft和ifft函数调用方式(1)Y＝fft(x)参数说明如果x是向量，则采用傅里叶变换来求解x的离散傅里叶变换；如果x是矩阵，则计算该矩阵每一列的离散傅立叶变换；(2)Y=fft(x，N)参数说明N是进行离散傅立叶变换的x的数据长度，可以通过对x进行补零或截取来获得。(3)Y=fft(x，[]，dim)或Y＝fft(x，N，dim)参数说明在参数dim指定的维上进行离散傅立叶变换；当x为矩阵时，dim用来指定变换的实施方向:dim=1，表明变换按列进行，dim=2表明变换按行进行。ifft的参数应用与fft完全相同。例4-1fft(x)函数的应用。X=[1423]Y=fft(x)运行结果：Y=10.0000-1.0000-1.0000i-4.0000-1.0000+1.0000i例4-2ifft(x)函数的应用。Y=[10.0000-1.0000-1.0000i-4.0000-1.0000+1.0000i]X=ifft(Y)运行结果：X=14232．fft2和ifft2函数调用方式(1)Y＝fft2(x)参数说明如果x是向量，则与一维快速傅里叶变换fft相同。如果x是矩阵，则计算该矩阵的二维快速傅里叶变换，数据二维傅里叶变换fft2(x)相当于fft(fft(x))，即先对x的列做一维傅里叶变换，然后对变换结果的行做傅里叶变换。例4-3fft2的应用。X=[2375;2465;2456;1234];Y=fft2(X);运行结果：Y=61.0000-14.0000+7.0000i-5.0000-14.0000-7.0000i0-7.0000i-3.0000+2.0000i4.0000-1.0000i-1.0000+2.0000i7.0000-2.0000+1.0000i1.0000-2.0000-1.0000i0+7.0000i-1.0000-2.0000i4.0000+1.0000i-3.0000-2.0000i例4-4ifft2的应用。对例4-3中的结果Y进行反变换。X=ifft2(Y)运行结果：X=23752465245612344.8基于DFT和FFT的频谱分析

4.8.1频谱分析的概念

频谱分析就是计算信号各个频率分量的幅值、相位和功率，以获得信号的频率结构以及各谐波和相位的信息。在实际工程中，大部分的信号都是模拟信号。对模拟信号x(t)进行频谱分析，就是计算信号的傅里叶变换。由于模拟信号及其傅里叶变换都是连续函数，不能用计算机实现，因此通常都需要通过时域采样将模拟信号x(t)转化为离散时间信号x(n)，然后利用DFT和FFT进行频谱分析。各频率分量幅值、相位以及功率的计算如下：4.8.2频谱分析的参数选择

在频谱分析中，参数的选择至关重要，这将影响频谱分析的准确性与精度。这些重要参数包括采样频率fs、频率分辨率Δf、最小的数据记录时间Tpmin和频域采样点数N。为保证采样后的信号不发生频率混叠，采样频率的选择必须遵循奈奎斯特采样定理，即在实际频谱分析中，采样频率的选择通常为信号最高频率的3~4倍。频率分辨率Δf频率分辨率（frequencyresolution）是频率域的采样间隔，即用DFT进行频谱分析时，能够分辨的两个频率分量之间的最小间隔。频率域的采样间隔为，因此频率分辨率Δf为：3.最小的数据记录时间Tpmin4.频域采样点数N在实际频谱分析中，频域采样点数就是进行DFT的点数，通常选择为2的整数次幂。例4-5对实信号进行频谱分析，信号最高频率为2.5kHz，要求频率分辨率小于10Hz，试确定以下参数：最短信号记录时间；最大采样周期；最少采样点数。例4-6已知信号x1(t)=sin2(30)t，x2(t)=sin2(120)t，x(t)=x1(t)+x2(t)，采样频率fs＝1000Hz，试利用MATLAB对信号x(t)进行频谱分析。解：MATLAB程序如下：t=0:0.001:0.6;x1=sin(2*pi*30*t);figure;subplot(3,1,1);plot(x1);

x2=sin(2*pi*120*t);subplot(3,1,2);plot(x2);

x=x1+x2;subplot(3,1,3);plot(x);

Y=fft(x);Pyy=Y.*conj(Y)/512;f=1000*(0:256)/512;figure;plot(f,Pyy(1:257))title('Frequencycontentofy')xlabel('frequency(Hz)')结果如图所示。例4-7假设信号含有3种频率成分，f1=20Hz，f2=25Hz，f3=35Hz，采样频率为100Hz，试用频谱分析的方法鉴别出每个频率成分。解：MATLAB程序如下：t=0:0.01:1;x=sin(2*pi*20*t)+sin(2*pi*25*t)+sin(2*pi*35*t)

;Subplot(2，1，1)

;plot(x)

;Y=fft(x，256);f=100*(0:(128))/256;Subplot(2，1，2)

;