《数字信号处理教程-MATLAB释义与实现》第四章

第四章

信号频谱的高效计算1信号频谱的高效计算4.1各种傅立叶变换及其关系4.2快速付立叶变换(FFT)4.3用FFT计算序列的频谱4.4时域采样中的频谱变换4.5连续信号的频谱计算4.6用反变换从频谱计算信号4.7用FFT计算能量4.8小结2（1）时域的周期性对应于频域的离散化。（2）

时域的离散化对应于频域的周期性。其主周期就是乃奎斯特频率范围[-π/T,π/T]。（3）

周期性的离散序列将对应于离散并周期性的频谱。即离散傅立叶级数(DFS)。离散有利于数值计算，但信号无穷延伸又不利于计算。（4）把时域和频域数据长度都限于主周期，并且使之相等，形成离散傅立叶变换(DFT)。它既离散，又长度有限，适合于计算机数值计算。

4.1各种傅立叶变换及其关系3各种傅立叶变换的特点变换名称时域信号（傅立叶反变换）频谱曲线（傅立叶变换）(连续)傅立叶变换(CFT)

连续信号连续频谱(连续)傅立叶级数(CFS)周期性，连续信号离散频谱离散时间傅立叶变换(DTFT)离散信号周期性，连续频谱离散傅立叶级数(DFS)周期性，离散信号周期性，离散频谱离散傅立叶变换(DFT)有限长离散信号(隐含周期)有限长离散频谱(隐含周期)4对离散傅立叶变换(DFT)，人们开发了可以高效地进行计算的方法，称为快速傅立叶变换(FFT)。人们想尽量利用FFT来解决其他类型信号的频谱计算问题。所以要充分弄清各种傅立叶变换之间的关系。本章就讨论这个主题。先把主要结果列出，其中有些结论已经讨论过，与采样定理有关的结论将在后面讨论，读者可先接受下来。目的是走通下图的路线，完成时频域傅立叶变换的数值计算。各种傅立叶变换及其相互关系5各种傅立叶变换及其相互关系

模拟信号时域采样周期延拓主值区间时域xa(t)——x(n)————x(n)·RN(n)|

|(FFT)频域Xa(jΩ)——X(jΩ)————X(k)·RN(k) CTFTDTFTDFS DFT ICTFT IDTFT 频域采样

IDFT6各种傅立叶变换及其相互关系(1).DFT与DFS的关系： 或 (4.1.7)(2).DFT与DTFT的关系

:采样插补

7各种傅立叶变换及其相互关系(3).DTFT与CTFT的关系:用采样定理或乃奎斯特定理建立关系。由CTFT求DTFT由DTFT求CTFT:频率泄漏可以忽略不计时近似解。

8各种傅立叶变换及其相互关系DFT与CTFT的相互关系:由CTFT求DFT只是取出其一部分，没有误差问题。而由DFT求CTFT要先经过由DFT求DTFT，再经过由DTFT求CTFT的两步。由离散求连续不一定有解，只有在信号的频谱足够窄，频率泄漏可以忽略不计的条件才可能近似由下式求得。

94.2快速付立叶变换(FFT)计算DFT的运算次数按N2快速增长。设N可以被2整除，把x(n)分成两个子序列x1(n)和x2(n)，则原序列的DFT可写成(设N1=N/2)：10快速付立叶变换(FFT)设它们的傅立叶变换分别为X1(m)和X2(m)，其周期为N1=N/2：则x(n)的傅立叶变换X(m)可表为：如果N=8,则N1=4,故X(m)周期是8,而X1(m)和X2(m)周期是4,即X1(4)=X1(0),X1(5)=X1(1),…。依次取m=0,1,…7，上式对应于右方的运算图。11快速付立叶变换(FFT)12快速付立叶变换(FFT)X1(m)和X2(m)是N1点的DFT，它们的计算又可以用类似方法化为两个更短的N2=N1/2点的DFT,一直分解下去，直到2为止，这就构成了上述FFT的全部运算图。粗算其中每一根线条代表一次乘法，线条的合成点代表一次加法。则每一级要N次乘法和加法。N=8时，需log28=3级，故共要24次乘法和加法。原来要N2=64次。若N=1024，需要10242次乘法，而用FFT,需分解为log21024=9级，只需1024×9次乘法，加快了100多倍。13快速付立叶变换(FFT)把上述运算次数的估计精确化：(1).每次分解的乘法次数为N，共log2N次， 乘法次数＝N·log2N(2).考虑到蝶形运算，又把乘法次数减半。在FFT运算图中，基本单元为如右图所示的蝶形结构，实际上不需要四次乘法，而只需要两次乘法即可完成。14蝶形运算节省一半乘法考虑到在算时，把m和N/2+m两点成组来进行，即构成上图的蝶形运算，就节省一半乘法15快速付立叶变换(FFT)用这些措施后总的乘法次数约为（当N很大时）：当N=1024时，结果为5120，与10242＝1048576相比，减小了近200倍。FFT除了提高速度，还要减少计算时所用的内存。理想的方法是实现原位计算，即每次运算的结果就放在输入数据的位置上，最后输出结果就放在原输入数据的位置上。这样所需的内存数目就是N个。为此，在具体实现FFT时，要遵循一些规则和技巧：16FFT如何节省内存(排序）(1).输入数据和输出数据的下标按二进制排序：时域抽取法(Decimate-in-time—DIT)—输入数据遵循倒序排列，则输出按顺序排列。如N=8,二进制顺序为[000,001,010,011,100,101,110,111]：即：0,1,2,3,4,5,6,7二进制倒序为[000,100,010,110,001,101,011,111]： 即：0,4,2,6,1,5,3,7故DIT输入数据按x(0),x(4),x(2),x(6),x(5),x(3),x(3),x(7)排列，输出下标即为顺序排列。17快速付立叶变换(FFT)(2).运算结构图按蝶形运算成对安排；每一段蝶形运算完成后，左边的一对数据就无用了，它们的位置就用来保存右边的一对数据。(3).运算从左到右，按段进行。每一段运算完成后，左一列数据就用不到了，全用右一列数据来取代。(4).各段全部运算完毕后，N个输入数据就换成了N个输出数据。其下标序号则由倒序变成了顺序。18快速付立叶变换(FFT)快速付立叶变换的其他变型和讨论：以上讨论的方法称为基2-DIT-FFT方法(1).仔细分析可以发现，N不必分解到2(基2)，基4的计算速度是最快的。(2).DIF-FFT方法与DIT-FFT形成对偶；(3).N不是2的幂次时可按素数分解,构成其他FFT算法.(4).FFT算法通常都用汇编语言编写，以进一步提高运算速度，所以各种软件系统都有FFT模块可供调用。一般不需自己编程。19快速付立叶变换(FFT)用MATLAB模仿FFT算法

functiony=myditfft(x)

m=nextpow2(x);N=2^m;%长度x对应的2的幂次m

%若x的长度不是2的幂，补零到2的整数幂

iflength(x)<N

x=[x,zeros(1,N-length(x))];end

%求1:N数列的倒序

nxb=dec2bin([1:N]-1,m);%求二进制顺序值

nxd=bin2dec(fliplr(nxd))+1;求十进制倒序值

y=x(nxd);%将x倒序排列作为y的初始值

（接下页）20快速付立叶变换(FFT)formm=1:m%将DFT作m次基2分解,

对每次作DFT运算

Nmr=2^mm;u=1;%旋转因子u初始化

WN=exp(-i*2*pi/Nmr);%基本DFT因子

forj=1:Nmr/2%本次间隔内的各次蝶形运算

fork=j:Nmr:N%本次蝶形运算跨越间隔 Nmr=2^mm;kp=k+Nmr/2;%蝶形运算的下标

t=y(kp)*u;%蝶形运算的乘积项

y(kp)=y(k)-t;%蝶形运算的加法项

y(k)=y(k)+t;%蝶形运算的减法项

end

u=u*WN;%修改旋转因子,

end,end21快速付立叶反变换(IFFT)快速付立叶反变换(IFFT):它和正变换是互成对偶的并有相似的形式:

两端取共轭：

可见x*(n)是X*(k)/N的傅立叶变换，故要求X(k)的反变换x(n),可先求X*(k)/N的FFT,再把求得的结果取共軛。22用MATLAB计算FFT和IFFTMATLAB中FFT函数的调用：正变换：X=fft(x);或 X=fft(x,N)N指定了FFT的长度。在默认条件下，它取x的长度。因为当N取2的幂次时，计算的速度最快，所以当x的长度不是2的幂次时，应尽量指定N。反变换用的是ifft函数，其调用方法与fft函数类似x=ifft(X);或 x=ifft(x,N)23数字信号处理器(DSP)中的FFTFFT进行的是重复性的乘积累加计算（MAC—Multiply-Accumulate-Calculation），普通PC机用的通用处理器需要用许多个指令周期才能完成一个MAC。所以开发了专门的DSP处理器，其主要特点是能高效地实现MAC运算。例如最新的DSP处理器可以在一个指令周期内完成一条MAC，并且在硬件中实现位倒序。在所有的DSP芯片开发系统中，FFT也都是标准的固件。

244.3用FFT计算序列的频谱

有限长序列的频谱计算：设已知位于n1~n2区间的x(n)，要求其频谱。在使用FFT时，必须用主值区间n=0,…,N-1的数据，FFT函数会产生N个数据，它们应该定位在下列频点上。

如果要求频谱关于零频率对称，可以利用X1=fftshift(X)函数。此时X1应分布在下列频点上。

25用FFT算序列频谱例4.3.1考虑长度为5的有限序列:x=[1,3,5,3,1],要求用FFT来计算其频谱。解：因为要求的是频率域中的连续频谱，故求出x的DFT和DTFT,进行比较。程序hc431rx=[1,3,5,3,1];nx=0:4;%给定原始数据N=length(x);D=2*pi/N;%求出及频率分辨率k=floor((-(N-1)/2):((N-1)/2));%求对称位置向量X=fftshift(fft(x,N));%求对称的FFT序列值subplot(1,2,1),plot(k*D,abs(X),'o:')%画幅特性subplot(1,2,2),plot(k*D,angle(X),'o:')%相特性图26用FFT计算序列的频谱曲线如下。注意plot语句具有线性内插作用，可近似得出连续频谱DTFT.27用FFT计算序列的频谱用FFT计算了频谱时，其频率样本数等于序列样本数，上例中于N=5，所以相邻两个频率样本点的间距为这个频率间距非常大，分辨率很差。改善分辨率的最好方法是增加数据的长度N。如果没有更多数据，也可以给输入序列补零。补了1019个零以后，得出连续光滑的频谱。不过改善的只是频谱的视在分辨率。

28用FFT计算序列的频谱例4.3.2考虑长度为11的矩形窗函数序列，计算其频谱。解：假如选N=20作为重复周期，则要在序列后面补9个零。使用FFT时，我们必须按N=20的周期延拓序列中取从n=0到19的主值部分，因此FFT的输入为

x=[ones(1,6),zeros(1,N-11),ones(1,5)]

在编写程序时，要准备给出不同的N进行比较。所以程序hc432的写法应能适用于不同的N。29序列频谱的程序hc432rC=[20,1024];forr=[1,2]N=C(r);D=2*pi/N;x=[ones(1,6),zeros(1,N-11),ones(1,5)];%主值区间的xk=floor(-N/2+0.5:N/2-0.5);%频率位置向量X=fftshift(fft(x,N));%求出x的对称位置FFTsubplot(1,2,r),plot(k*D,X)%画实频特性图end30用FFT计算序列的频谱由图，N=20给出了频谱的大体形状，但分辨率差，N=1024的结果则接近于精确的结果。

31用FFT计算序列的频谱长N的有限序列的频谱可以用N点FFT方便地求出。存在的问题是：频谱是定义在奈氏频率范围内所有的ω上的，而FFT只计算样本点上的值。因此，要由插值才能求出全部点的频谱。较简单的插值方法是在输入数据的尾部补零，以改善频率分辨率。如果一个有限序列位置向量不在主值区间，比如在某些 ，则这些负时间的部分必须移到补零后的尾部，使位置向量位于[0:N-1],然后再用FFT。否则得到的相频特性将是错误的。

32用FFT计算无限序列的频谱当把一个无限序列截断成为有限序列，会导致泄漏和波动，算出的频谱与无限序列的频谱有误差。因此要确定截断的长度至少应有多长，才能保证其频谱足够准确。

但原无限序列的频谱又是未知的，怎么知道误差？所以就采用另一种思路：不断地倍增截断的长度，把相邻两次计算的结果进行比较。如果其间的误差已远小于算出的频谱峰值，则认为该结果已经接近于（收敛于）精确值。33用FFT计算无限序列的频谱例4.3.3计算下列序列x(n)的频谱,精度要求为频谱峰值的1%。解：这是理想题，它的频谱已经解析地算出为：现在要问，应该截取多长的序列，才能使计算出的频谱达到给定的精度。以后再考虑不知道频谱解析解的情况。34无限序列频谱的算例4.3.3rN=32;dw=2*pi/N;n=0:N-1;x=0.7.^n;k=n;X=fft(x); %求出32点数据的DFT数值解%解析法得出的32点频谱XtXt=1.0./(1-0.7.*exp(-j*k*dw));%求数值解和解析解的最大误差e=max(abs(abs(X)-abs(Xt)));Xm=max(abs(X));%求序列的最大幅度pe=e/Xm*100 %求出最大相对误差（百分数）

35无限序列的频谱计算程序算出的百分数误差pe是0.0011,相当于十万分之一。可见取32个样本足够了。现在设不知道准确频谱。要选定最小的N值，使相对误差不超过β%。这时，要逐次倍增N值，比较相邻两次计算结果的误差，以这个误差作为判断的标准。使得用N=2a个数据点计算出的频谱，与用N/2个数据点的计算结果的误差小于峰值的β%，比较频谱的误差必须在相同的频点上进行。因为点数N变了，同样的频点在相邻两次计算中的下标不同。

36比较N1=N和N2=2N的计算频点。用N1点算出的频谱X1位于下列频点上

(4.3.4)而用N2点算出的频谱X2，则位于下列频点上

(4.3.5)令 ，解出k2=2k1。必须在的这些频点上比较X1和X2的幅度。按照这个原理来编写程序hc424。

无限序列的频谱计算37计算程序例hc434ara=1;b=100;beta=1;%给定初始数据whileb>beta%判断是否应结束循环运算

N1=2^a;n1=0:N1-1;%确定数据长度N1x1=0.7.^n1;X1=fft(x1);%长度N1的x1及其FFT1N2=2*N1;n2=0:N2-1;%数据长度加倍

x2=0.7.^n2;X2=fft(x2);%长度N2的x2及其FFTk1p=0:N1/2-1;k2p=2*k1p;%两序列对应k2=2k1d=max(abs(X1(k1p+1)-X2(k2p+1)));%对应点的误差

mm=max(abs(X1(k1p+1)));%X1的最大幅值

b=d/mm*100;a=a+1;%求相对误差,长度加倍end38在此程序中，如果设β=1，程序将得到N2=32，b=0.3323，满足了b<β的要求。这时的b基本上反映的是N=16时的误差，也就是说取N=16就够了。不过既然已经按N=32算出了较精确的结果，也没必要再退回去取N=16重算了。比较相邻两个N值的计算误差列表如下：

程序hc434的计算结果N1/N24/88/1616/3232/6464/128

误差b%24.015.76480.33230.00111.22e-8394.4

时域采样中的频谱变换

时域采样定理(Nyquist定理)需要弄清经采样后的序列与原始的模拟信号在时域和频谱方面的数学关系。先把采样的过程看作模拟信号通过一个电子开关S；再考虑增益而改为乘法器。见图4.1.1(a)和(b)。把开关信号等价成一宽度为τ，周期为T、面积为1（因而幅度为1/τ）的矩形脉冲串 ，采样信号就是xa(t)与 相乘的结果。如果让电子开关合上时间τ→0，则形成理想采样，此时上面的脉冲串变成单位冲激串40连续信号采样中的频谱变换采样过程用电子开关模型(a)和乘法器模型(b)来表示。后者将采样过程看作把xa(t)与宽度为τ,周期为T,面积为1，因而幅度为1/τ的矩形脉冲串PT(t)相乘，当宽度

τ趋向于零时，这些脉冲趋向于δ函数。41连续信号采样中的频谱变换因此连续信号与其采样序列间的时域关系表为：对等式两端求傅立叶变换，由于 而即单个脉冲具有无限宽的均匀频谱。当单个脉冲构成等间隔脉冲序列时，42连续信号采样中的频谱变换它的频谱是一个等幅度、周期性的脉冲频谱。因为时域相乘对应于频域卷积：δ函数的一个特性是：它与任何函数的卷积就等于该函数。一串间隔为Ωs的δ函数与X(jΩ)卷积的结果是43连续信号采样中的频谱变换这是采样序列频谱和原模拟信号频谱的关系，注意两者也差了一个时间量纲T，与时域相同。上式表明，采样信号的频谱是一串原模拟信号的频谱无限次平移的叠加。平移的间隔为采样角频率Ωs=2πFs。可以看出，知道Xa(jΩ)求X(jΩ)是有唯一解的，而知道X(jΩ)求Xa(jΩ)则不是唯一的。44连续信号采样中的频谱变换设模拟信号xa(t)是带限信号，最高截止频率为Ωc，其频谱Xa(jΩ)如

(a)图。Pδ(t)的频谱Pδ(jΩ)如图(b)。是一连串间隔为Ωs的脉冲频谱串。x(t)的频谱X(jΩ)是(a)和(b)的卷积。Ωc<Ωs/2时如图

(c)，原模拟频谱互相分开。Ωc>Ωs/2时如图

(d)，原模拟频谱互相交叠。45连续信号采样中的频谱变换设模拟信号xa(t)是带限信号，最高截止频率为Ωc，其频谱如图4.4.2(a)所示。Pδ(t)的频谱如图4.4.2(b)所示。那么按照（4.4.7）式，x(t)的频谱如图4.4.2(c)所示。如果基带频谱较窄，满足Ωc<Ωs/2，基带谱与它的周期延拓的谱没有重叠，此时的序列频谱就如图4.4.2(c)所示的周期频谱。如果让这个采样信号经过理想低通滤波器G(jΩ)，就能得到输出y(t)，不失真地提取原模拟信号，如图4.4.3所示。此时，(4.4.7)式成为：46连续信号采样中的频谱变换对离散序列进行低通滤波[如图

(a)]复原连续信号的可行性：如果原模拟信号的频谱比采样频率之半低，即Ωc<Ωs/2，如图

(b);取低通滤波器的带宽等于乃奎斯特频率，如图

(c);则滤波输出的频谱就是原模拟信号的频谱如图

(d);47连续信号采样中的频谱变换(4.4.8)式说明，在没有混叠的情况下，序列的频谱等于原始模拟信号的频谱除以采样周期。如果选择采样频率Ωs低，或者说信号最高截止频率Ωc高，以至Ωc>Ωs/2，或者fc>Fs/2，则基带谱与把它作Ωs的周期延拓形成的谱就会发生重叠，称为频谱混叠现象，如图4.4.2(d)所示。这种情况下，再用图4.4.3所示的理想低通滤波器对它进行滤波，得到的将是失真了的模拟信号。

48连续信号采样中的频谱变换采样定理(Nyquist定理)：

采样信号频谱是原连续信号的频谱以采样频率Ωs为周期进行周期延拓并叠加形成的，用公式（4.4.7）表示。

Ωs/2称为折叠频率或乃奎斯特频率。如信号带宽大于折叠频率，Ωc>Ωs/2，则采样信号x(t)的频谱会是原模拟信号频谱折叠相加，称为混叠或泄漏，造成频谱的失真。实际中需根据模拟信号的最高频率Ωc，按Ωs>2Ωc的要求选择采样频率Ωs(Fs=Ωs/2π

）。49模拟信号的重构重构模拟信号是采样的逆过程——从离散到连续。设序列x(n)的频谱为X(jω)，根据公式(4.4.7)，它与模拟频谱Xa(jΩ)关系如下：如果信号带宽小于折叠频率，则可以用带宽为[-π/T,π/T]的理想滤波器完全无失真地恢复原来的低频模拟频谱。在时域也可完全无失真地恢复原来的低频模拟信号波形。50模拟信号的重构由离散序列重构模拟信号的过程可以分解成两步。(1).首先经过一个脉冲变换器，将离散序列乘以T，使它变成一个等价的模拟信号；(2).然后再让它通过一个低通滤波器。滤波器的输出就是模拟信号。下面将分理想方法和实际方法两种情况讨论。

51模拟信号的重构

理想方法：

(1).脉冲变换器把每个数字样本x(n)变成相应权值的δ函数，即产生如下的时间函数xs(t)(2).然后xs(t)让通过一个幅频特性如图4.4.3(c)的理想重构滤波器。其数学形式为：

52模拟信号的重构对上述理想频率特性作傅立叶反变换（例3.2.8),得到它的脉冲响应：右图是它的波形，下图是使它右移一拍。可以看出，它是非因果的，所以是不可实现的理想滤波器。53模拟信号的重构滤波后的输出ya(t)应为xs(t)和h(t)的卷积：它是理想滤波器的脉冲响应平移加权的叠加，如右图所示。54模拟信号的重构由于理想脉冲响应h(t-nT)在t=nT点都等于零，叠加的结果对除中心样本点外其它样本点的值没有影响。所以输出的模拟信号ya(t)将经过所有的样本点x(n)，它就是样本序列的包络。理想方法的两个步骤都是实际上无法实现的。第一步中的脉冲函数δ(t)在工程中无法做到；第二步中的理想滤波器的脉冲响应又是非因果过程。所以理想方法只有理论意义，它给出了一种理论计算的方法和努力追求的目标。55模拟信号的重构实际方法：(1).用零阶保持器（zero-order-holder—ZOH）代替脉冲变换器。它把序列保持一个采样周期T，直到下一个样本的到来。（图1.1.1)这样，就相当于把信号乘了T，并使之成为模拟信号；

(2).零阶保持器自身也有低通滤波作用，因为它的脉冲响应为：56模拟信号的重构用傅立叶变换求得并画出零阶保持器的频率特性：57模拟信号的重构零阶保持器的幅频特性的通带约为π，即乃奎斯特频率。这是很好的特性。从时间波形而言，它大体相当于图2.1.1(b)，但保持波形右移了半拍。所以它基本反映了样本序列的包络，不过在时间上延迟了半个采样周期。这也可以从它的相频特性上看出。作为滤波器，零阶保持器的过渡带太缓，高频区又有较大的波动。反映在波形中有高频的尖峰，通常还不能满足工程的指标，往往还要加上后滤波器。58预滤波作用的定量分析按两种情况作对比分析：(1)在图4.4.7(a)中，连续信号xa(t)经过预滤波后为xp(t)，然后把它采样后变成x1(nT)。再用一个带理想低通滤波器的D/A变换器把它转换回连续信号xa1(t)。(2)在图4.4.6(b)中，连续信号xa(t)不经过预滤波而被采样，得到的x2(nT)也立即再用同样的D/A变换器把它转换成连续信号xa2(t)。比较误差xa1(t)-xa(t)和xa2(t)-xa(t)

哪个大？59预滤波作用的定量分析60预滤波作用的定量分析用误差平方积分E1和E2和帕塞瓦尔定理：由于x1(n)没有混叠，在奈奎斯特主频区间有经过理想重构后，只保留主频区间频谱不变，别处为零61预滤波作用的定量分析平方积分E1中[-π,π]的一段积分为零：可写成对于E2,，由于存在着频率混叠，在|Ω|≤π/T的范围内，Xa2(Ω)将不等于Xa(Ω)，因此平方积分E2可分成三段来写：62预滤波作用的定量分析与E1相比，其首尾两段相同，所以可写成这证实了信号经过预滤波后再作数字处理的比不经过预滤波的误差要小。因为这个误差主要取决于采样后造成的频率混叠。预滤波器的带宽应当取得比乃奎斯特频率π/T小。634.5连续信号的频谱计算对于实际中遇到的大多数连续信号，都不能采用解析方法求频谱。而必须采用计算机辅助的方法。第一步是采样，把连续信号离散化；第二步是截断，取有限的数据。不然就用不成计算机。所以，不仅是实时处理信号要用数字信号处理的理论，非实时地分析研究大部分连续信号照样要用数字信号处理的理论和方法。64设xa(t)为绝对可积的连续时间信号，x(nT)是它的采样序列，如果T足够小，有数字频谱X(jΩT)=X(jω)可以用FFT求得，所以就可由上式求得连续时间信号的频谱。计算中的关键问题是如何选择采样周期T、数据长度N和/或信号持续时间L。三个变量之间以L=TN相关联，只有两个独立可选。连续信号的频谱计算65连续信号的频谱计算在参数选择时，要考虑如下三个问题：1．频率混叠：如果不是有限带宽的，则采样周期T的选择的主要依据是，取足够小的T，使得由时间采样所造成的频率混叠小到可以忽略。2．频率分辨率：连续信号频谱Xa(Ω)是定义在所有Ω上的，当使用FFT来计算时，只计算出了它在N个数字频率上的样本值X(jω)，ω=2kπ/N,k=0,1,…,N-1。数字频率分辨率为dω=2π/N，选择大的N可以改善数字频率分辨率，获得频谱细节的更多信息。663．截断效应：如果信号是无限长的，就必须把它截断到长度L=NT，截断相当于加窗，它会引起吉布斯效应（波动），也会把窗函数的频谱引入信号频谱，造成混叠，需要考虑其误差的问题。此外，模拟频率分辨率为D=dω/T=2π/L，它与L成反比，选择大的信号持续时间L=NT可以改善模拟频率分辨率，得知频谱的细节。

连续信号的频谱计算67连续信号的频谱计算由于上述三个问题，当用FFT来计算的频谱时必须小心。很明显，采样周期T应该选得小些以减少混叠，而N要选得足够大来提高分辨率；如果N是确定的，为了减小频率混叠就要减小T；为了减小截断效应，需要增加L=NT，所以要增大T；这是两个互相矛盾的要求。此外，取很小的T和很大的N往往又是不必要的。因此建议在计算与频谱时，采用以下的步骤来选择T和N。68连续信号的频谱计算1．先初步选择时间纪录长度L，使得0到L之间包括了大部分非零的xa(t)，然后用逐次减小T和加大N的步骤来选择周期T，使得时间采样造成的频率混叠可以忽略不计。2．在选定T以后，进一步选择L。即把L(和N)增加一倍。并与用原来的结果进行比较。因为两个频谱所用的T相同，因此，两个频谱之间的差别是由于截断效应不同（即L不同）引起的。如果两者的误差不大，这个L就可以接受；反之，再把L加倍，直到用两个相继的L算出的频谱非常接近为止。69连续信号的频谱计算例4.5.1考虑连续时间信号:用FFT计算其频谱。解：先选tmax=L=10, ，所以L=10并没有完全包括信号的主要部分。但时间区间[0,10]还是能大体反映t>0的信号的主要频率分量。再随意地选择N=5，得到T=L/N=2。然后编写计算频谱的程序hc451,核心语句为：70连续信号频谱计算例4.5.1T=？;N=？;D=2*pi/(N*T);n=0:N-1;x=exp(-0.1*n*T); ％求出离散序列X=fftshift(fft(x));％用求对称离散频谱，Xa=T*X； ％乘以T求连续频谱k=floor(-(N-1)/2:((N-1)/2);％位置向量也取对称Xa(1) %求乃奎斯特频率处的频谱值plot(k*D,abs(Xa)) %绘图

程序中有几个地方需加说明：71连续信号频谱计算例4.5.1(1).从模拟信号求连续频谱，经过的步骤是：连续信号xa(t)→(采样)离散序列x(n)=xa(nT)→(FFT)序列频谱X(k)=fft(x)→(乘T)连续信号频谱Xa(kD)=TX(k)，这反映在2～4行程序中.(2).检验混叠，用的方法是在折叠频率上的频谱幅度足够小。以判断参数T选择是否合理。故求对称频谱向量的第一项X(1)的值。(3).选择N,用不同T和N算得到的幅频特性画在图中。注意六张子图上的频率坐标是相同的，但它们的乃奎斯特频率边界[-π/T,π/T]不同。72连续信号的频谱计算73连续信号的频谱计算T的选择：可以用乃奎斯特频率处的幅特性来评价混叠的严重程度。图4.5.1(a)中T=2,乃奎斯特频率边界是

[-1.57,1.57];，图(b)的T=1,边界是[-3.14,3.14];（超出了图中横轴的边界）。可以看出，在图4.5.1(b)的边界(±3)处，有一个相当大的非零值Xa(1)=0.3319。说明T=1时，频率混叠仍不能忽略。继续多次把T减半和把N加倍（保持L不变）进行计算，直到T=0.1,N=100，这时的乃奎斯特频率已经是±31.4，其上的Xa(1)值已减小到0.0318，峰值是6，相对误差小于0.032/6=0.4%。这说明，选T=0.1，频率混叠可以忽略不计。

744.5连续信号的频谱计算N的选择：图4.5.1(e)用实线画出了T=0.1,N=200的计算频谱，与图4.5.1(d)中T=0.1,N=100的结果进行比较。这里只画出了[-3,3]的频率范围，可见它们之间的区别很大，说明截断效应相当明显，有必要采用更大的L。图4.5.1(f)给出了用T=0.1,N=400的计算频谱。与图4.5.1(e)相比仍有相当差别。如果再画出T=0.01,N=800的结果，它们的误差在[-3，3]频段上就很难分辨了。因此取L=0.1×800=80秒已经大得足以避免截断效应的影响。

75截断长度和补零的影响Xa是t=0到L=40，采样周期为0.1秒所得的全部400个数据，即xa=exp(-0.1*(1:N)*T);xa1是只取xa的前N2=100个数据,而后面N-N2个数据补以零的数组xa1=[exp(-0.1*[1:N2]*T),zeros(1,N-N2)];xa和xa1的波形及其频谱见图4.5.2。由于截断，频谱Xa1有波动。即窗函数的吉布斯效应。解的结果为：可以从中看出波动分量。76截断长度和补零的影响77截断长度和补零的影响可见，序列补零可以得到高的分辨率，经FFT得到很密的频谱，但频谱的值不见得精确，因为它没有增加信息。在实用中，如果可以得到的更多数据时，应该补充数据而不是补零。补零可能反而给出了频谱中不正确的细节。在本例中，时域截断是计算中人为地引入的，其频域波动也就是人为引入的。所以首先要解决的是原始时域数据改善的问题，而不是去补零来提高分辨率。

78连续信号的频谱计算为了用一个程序而，用不同的参数画出图4.5.1的多幅子图。可以用for循环如下：T0=[2,1,0.5,0.1,0.1,0.1];%各次的T,编成向量L0=[10,10,10,10,20,40];%各次的N,编成向量forr=1:6%循环计算六次

T=T0(r);N=L0(r)/T0(r);%顺序调用T及N（下接程序hc451核心语句）subplot(3,2,r),%在第r个子图位置plot(k*D,abs(Xa))%绘图end79连续信号的频谱计算例4.5.2用FFT计算下列连续信号的频谱。解：本例的计算原理同上，因为信号由两个频率很接近的正弦信号组成，所以重点是计算时保证频谱的分辨率D,它取决于截断信号的长度L,D=2π/L。计算的程序为hc452,80连续信号的频谱计算T0=[0.6,0.15,0.15,0.15];N0=[256,256,256,2048];forr=1:4(类似于hc451的计算连续信号频谱核心语句）subplot(2,2,r),plot(k*D,abs(Xa))%画幅频特性end上述四组T0和N0对应于L取[153.6,38.4,307.2]三种数据.其频率分辨率为D=[0.0409,0.1636,0.0205]。由于两个正弦频率之差为0.1，所以只有T0=0.15和N=2048能同时保证较小的混叠和较高的分辨率。8182连续信号的频谱算例例4.5.3计算下列信号的频谱设a=5，这个信号是有限时间信号，它是实的和偶的，因此它的频谱也一定是实的和偶的。先任选T=1，则在|t|≤5的h(nT)中总共有11个样本点。因此，把时间样本数选成奇数，并且要按循环对称的概念排列在0到N-1的主值下标区间内。编出程序hc45383连续信号的频谱计算T=1;M=5/T;N=2*M+1;n=1:M;%原始数据D=2*pi/(N*T);%频率分辨率hp=sin(2*n*T)./(pi*n*T);%生成正序列h=[2/pihpfliplr(hp)];%构成循环对称序列H=T*fftshift(fft(h));H(1)%求xa的对称FFT，k=floor(-(N-1)/2:(N-1)/2);%频率下标向量

subplot(1,2,1),plot(k*D,H)%绘图

程序运行结果说明，N取11时分辨率太低，取T=0.5,用1024点FFT,(相当于后面补零到1024）可得到较好结果。84连续信号的频谱计算85连续信号的频谱计算周期信号的频谱计算：从理论上说，如果一个连续时间信号是周期性的，它必定有无限的长度，计算出的频谱将不会收敛。然而，用计算机计算时，序列总是有限长的。因为我们遇到的都是经过截断了的周期信号。从算法上看，完全可以用上节中的方法和程序，只是在计算结果上，频谱会有向有限个频点上集中，并成为脉冲的趋势。

86连续信号的频谱计算例4.5.4计算定义在全部t上的xa(t)=cos5t的频谱，它的理论频谱是 ,它包含了权重为π的位于Ω=±5上的两个脉冲函数。

解：在计算机计算中，正余弦函数必须截断为有限长度L。信号＝cos5t的带宽限制于5。纯理论地看，只要采样周期小于π/5=0.63秒，就不会发生频率混叠。然而如果把cos5t截断为长L的信号，则它的频谱就不再是有限带宽了，所以必须采用更小的采样周期，任选T＝0.1，并选N＝50，得到L＝TN＝5。按此来截断信号。

87连续信号的频谱计算程序hc454。N=input('N=');T=0.1;n=1:N;%原始数据D=2*pi/(N*T);%频率分辨率xa=sin(5*n*T);%生成有限长的正弦序列Xa=T*fftshift(fft(xa));Xa(1)%求xa的对称FFT，k=floor(-(N-1)/2:(N-1)/2);%频率下标向量

plot(k*D,abs(Xa)) %绘图

程序运行的结果见下图。88连续信号的频谱计算89连续信号的频谱计算把N作为用户自选的参数。选N=50,100,500及628所得的计算结果在图中画出，在±5处都出现了两个尖峰。结论是：如果用FFT计算出的幅频谱，包含了窄的尖峰，而N每加一倍，它也大体上加倍，这就说明这个连续时间信号包含一个周期分量。因为尖峰窄了以后，很难准确采样在峰值处，所以图上显示的峰值，并非真正的最大值。从理论上说，只有长度L恰好是周期的整数倍时，可以得到准确的峰值。

90连续信号的频谱计算例4.5.5求下列周期信号的频谱xa(t)=-1+2sin0.2πt-3cosπt解：它的最高频率是Ωm=π弧度／秒，因此采样周期T必须小于π／Ωm＝1。任选T＝0.2，在图（a）中画出N＝50时的幅频特性，而在图(b)中用N＝1000。

从图(b)中也可大致看出，在ω＝0，±0.628和±3.14的三个高度之比为2：2：3，由于信号的直流分量是它的零阶傅立叶级数的一半，可见原信号的三个频率分量幅度之比为1：2：3，91连续信号的频谱计算92连续信号的频谱计算用FFT从连续信号计算其频谱的步骤概括如下：连续信号xa(t)―→(采样)

离散序列x=xa(nT)―→(FFT)

序列频谱X(k)=fft(x)―→(乘T)信号频谱Xa=T‧X1.采样前的xa本来是解析式，但在MATLAB中和x同样用数组表示，所以在程序中用同一变量表示；2.序列样本经FFT后变为序列频谱，其频率混叠误差由折叠频率处的幅特性决定。该处的幅特性愈小愈好。3.控制混叠误差和频率分辨率的参数是T,N和L。4.序列频谱必须乘T才等于信号频谱。这也表明信号和它的采样序列物理上不等价。93循环计算中对应频点的确定

当要靠计算机自动选择参数N,T或L时，需要比较相邻两次计算中对应频点上幅特性的误差。本节讨论在N,T或L发生变换的条件下，如何找到对应的频点Ωk？

因为当相邻两次计算中，N、T或L发生变化，必须要在同样的Ωk下进行比较。94循环计算中对应频点的确定令得到：只有满足(4.5.7)式的整数k1和k2，才是可以进行频谱比较的对应点。这是从数学上说，如果要得到更形象的结果，可以用由程序hc456生成的图4.5.7来观察。

95循环计算中对应频点的确定96循环计算中对应频点的确定以子图一作为标准，T=0.1秒，N=16，信号长度L=1.6秒，乃奎斯特频率范围为 ；内分16个点。把T减小一半，N不变，L也减小一半，得到子图二。其乃奎斯特频率范围扩大了一倍，成为 内分16个点，间隔大了一倍。故只有8个点能与标准情况对应。把T减小一半，同时把N加大一倍，因此L不变，得到子图三。其乃奎斯特频率范围也扩大了一倍，分成32个点，故中间16个能与标准情况对应。T不变，N加大一倍，L也加大一倍，得到子图四。其乃奎斯特频率范围不变，分成32个点，间隔D缩小一半，每隔一个点有一个点与标准情况对应。974.6用反变换从频谱计算信号

从连续频谱计算信号的步骤应为求频谱过程的逆：信号频谱Xa(Ω)→(频率域采样除T)

离散序列频谱X(k)=Xa/T→(IFFT)时间序列x(n)=fft(X)→(取包络)连续信号xa信号频谱必须采样并除T才等于序列频谱；序列频谱经IFFT后变为时间序列，其时间混叠误差由时间序列中幅度最小的区域决定。该区的幅度愈小愈好。3.控制混叠误差和时间分辨率的参数是T,N和L。4.

xa是x的包络，可以用plot(nT,x)画出。但作为变量，xa在MATLAB中和x用同一变量表示。98频率域采样定理

X(jω)

的IDTFT是有限长序列x(n)；其频率域采样序列X(k)(以后记作Y(k))的IDFT是无限长周期序列的主值y(n)。计算机不能做IDTFT，只能做IDFT，因而只能得到y(n)，得不到x(n)。从概念上说，X(jω)中包含的信息是完全的，Y(k)只是其部分样本，是不完全的。所以，知道X可以确定Y；知道Y未必能完全确定X。对应于时域信号。知道x(n)可以确定y(n)；知道y(n)未必能完全确定x(n)。频率域采样定理从数学上求出两者的数学关系。99证明的基本思路是：写出x(n)=IDTFT[X(jω)]；写出y(n)=IDFT[Y(k)]；通过Y(k)=X(jωk)=X(j2πk/N)把两者关联起来；频率域采样定理的证明100注意 ,即周期为N,把x(n)分成 长N的段。段号为r，段内下标为n1Є[0,N-1]，则x(n)=x(n1+rN)，r=0,±1,±2,…

（4.6.4）中n1相同的点的乘子 相同，可以合并得此式与(4.6.2)式恰成傅立叶变换对，故有可见，y(n)是无限个x(n)平移相加的结果。101由（4.6.7)，如果已知x(n)，可以将它无限次平移N相加而得到y(n)。如果已知y(n)，则未必可能求出x(n)。只有当x(n)在[0,N-1)主值区间外取值为零，即此时x(n)平移后不互相重叠,即无时间泄漏。

x(n)=y(n)(n=0,1,…,N-1)这就要把N取得足够大，即频率域采样足够密。频率域采样定理102从频谱计算离散时间序列为了利用IFFT函数，输入的频率样本点必须合乎IFFT的要求。(1).IFFT输入的N个样本应是数字频率ω＝ΩT,它的最大值不得超过2π。假如频谱的最大频率为Ωm，采样周期T必须满足：式中，Fs为采样频率(赫兹)。T或Fs确定后，才能把给定的频谱X（Ω）映射到数字频率坐标上来，成为X（ω）。频率分辨率为

Δω=2π/N或D=Δω/T=2π/(NT)103从频谱计算离散时间序列(2)通常序列的DTFT是用对称于原点的正负频率表示，频率下标取为

k＝floor(-(N-1)/2:(N-1)/2)

而IFFT的输入必须从k=0到N-1，因此，频率样本点中k<0的部分（用kn表示）必须平移到k≥0部分（用kp表示）的右面。

(3)若待求的时间序列是实的，则输入频谱X(Ω)应具有共軛对称性。X(k)应具有循环共軛对称性。所以把k作上述移动不影响其对称性。

104从频谱计算离散时间序列由上节，频率采样会引入时间混叠。如果时间混叠不能忽略，那样算出的结果是无用的。检查时间混叠的方法如下：如果算出x1(n)在所有的n处都明显地不等于零，说明时间混叠相当大，必须另选大些的N。当算出的x1(n)在n的某个范围内实际上等于零了，即IDFT存在着(4.6.8)式中含零的部分，就认可这个N。105从频谱计算离散时间序列例4.6.1设某时间序列的频谱为要求用IDFT数值计算求出它的时间序列。解:因为题中用的是数字频率而没有涉及T，这就是默认T=1.先取N=3，输入频率数组必须是把2π等分为3份的数字频率。然后再取N=10，用同一程序，对计算结果进行比较。106从频谱计算离散时间序列

程序hc461T=1;N=3;D=2*pi/N; %N及频率分辨率kn=floor(-(N-1)/2:-1/2);%负频率下标向量kp=floor(0:(N-1)/2); %正频率下标向量w=[kp,kn]*D;%形成主值区频率排序%按新的频率排序输入数字频谱X=2-exp(-j*w)+exp(-j*2*w)+exp(-j*3*w);x=ifft(X);%对X求IFFT,得出循环对称的xstem(T*[0:N-1],x)%绘图107从频谱计算离散时间序列108从频谱计算离散时间序列当N=3时，因为N小于时间序列长度，产生了时间混叠，其标志是序列中没有取零的点或点列，结果不正确；当N=10>4时，结果是正确的，其标志则是在算出的IDFT序列x(n)中，包含一段如(4.6.8)式所指的全零的部分，

109从频谱计算离散时间序列例4.6.2给出在频率范围[-62.8,62.8]弧度/秒间的频谱。求出它的时间序列。解：本题给出的是模拟频率Ω，Ωm=62.8取T≤π/Ωm=0.05秒。即采样频率为20赫兹或40π。N取8和20两种情况。

110从频谱计算离散时间序列wm=62.83;T=pi/wm; %选择采样周期TN0=[8,20]; %设定两种N值forr=1:2%作两次循环计算

N=N0(r);D=2*pi/(N*T);%取N,求频率分辨率

k=[0:N-1]+eps;%求频率下标，偏移微量

X=sin(0.275*k*D)./sin(0.025*k*D);算出频谱

x=ifft(X);%求X的IDFT，得出x(n)subplot(1,2,r),stem([0:N-1]*T,x)%绘图end111从频谱计算离散时间序列112从频谱计算离散时间序列N=8时得到左图，其中的8个数据中没有一个为零，可见存在着频率采样所造成的时间混叠，再取N=20进行计算。结果画在右图中，算得的序列中都有1