第5章连续时间信号与系统的频域分析






对于连续时间信号的频域分析，本章将首先介绍周期信号傅里叶级数的数值计算方法，再介绍非周期信号傅里叶变换的符号函数求解以及数值计算，最后验证部分傅里叶变换性质。
对连续时间系统进行频域分析的前提是该系统是线性时不变的，只有这样零状态响应yzs(t)才等于输入信号x(t)卷积单位冲激响应h(t)，从而零状态响应的傅里叶变换Yzs(jω)才等于输入信号的傅里叶变换X(jω)乘以单位冲激响应的傅里叶变换H(jω)，用H(jω)对系统进行频域分析才有意义。与h(t)一样，H(jω)仅取决于系统本身，与输入无关，是表征系统特性的一个重要物理量。





5.1周期信号的傅里叶级数
5.1.1傅里叶级数的计算


工程中主要采用指数形式的傅里叶级数对周期信号进行频域分析。周期为T的周期信号x(t)，若满足狄利克雷条件，在任意［t0,t0+T］区间，可用虚指数信号集ejkω0t,ω0=2πT,k=0,±1,±2,…精确分解为以下形式的傅里叶级数



x(t)=∑∞k=-∞Xkejkω0t
（5.1）


其中，


Xk=1T∫t0+Tt0x(t)e-jkω0tdt
（5.2）



在已知周期信号x(t)的数学表达式且能计算出积分的条件下，可以用式（5.2）精确计算傅里叶系数。但当数学表达式非常复杂，无法得出积分结果，或者写不出周期信号的数学表达式时，可以先对周期信号在一个周期内进行抽样，再进行数值计算。
下面讨论利用数值计算得到连续时间周期信号傅里叶级数的方法。设一个周期内刚好抽样得到了N个点，则抽样间隔Ts=TN，抽样得到的序列可以记作x(t0+nTs)，式（5.2）可表示为


Xk≈TsT∑N-1n=0x(t0+nTs)e-jkω0(t0+nTs)
（5.3）



例5.1.1计算如图5.1.1所示的周期矩形脉冲信号第-10~10项的指数形式傅里叶系数。


图5.1.1连续时间周期矩形脉冲信号


解： 可通过for循环直接得到Xk，MATLAB源代码如下

close all; clc;clear all;

T = 4;%周期 

N = 400;%一个周期抽样点数

dt = T/N;%抽样间隔

t=-T/2:dt:T/2-dt;%时域自变量

x=2*rectpuls(t,2);%x(t)的赋值

w0 = 2*pi/T;%基频

k = -10:10;%需要计算的Xk的项数

Xk = zeros(size(k));

for m = 1:length(k)

for n=1:N %MATLAB中数组从第1项开始

Xk(m) = Xk(m) + dt/T*x(n)*exp(-j*k(m)*w0*t(n));%按公式计算Xk

end

end

stem(k*w0,Xk,'filled'); xlabel('＼omega'); title('Xk')

因为Xk为实数，所以在一幅图里画出了频谱，结果如图5.1.2所示。将该结果与Xk的理论值Sakπ2进行对比，结果如果图5.1.3所示，增加抽样点数，数值解与理论值将会越来越相近。


图5.1.2连续时间周期矩形脉冲信号的频谱




图5.1.3Xk的数值解与理论值之差


因为MATLAB是基于矩阵的运算，前面用for循环计算Xk的方法虽然原理直观，编程简便，但耗时较多。可以将for循环改为矩阵向量相乘，降低计算复杂度。先将式（5.3）中Xk表示成向量相乘



Xk≈TsTe-jkω0t0e-jkω0(t0+Ts)…e-jkω0(t0+(N-1)Ts)
x(t0)
x(t0+Ts)
︙
x(t0+(N-1)Ts)


（5.4）



进一步地，将Xk表示成向量形式，从而



Xk1
Xk1+1
︙
Xk2≈TsTe-jk1ω0t0e-jk1ω0(t0+Ts)…e-jk1ω0(t0+(N-1)Ts)
e-j(k1+1)ω0t0e-j(k1+1)ω0(t0+Ts)…e-j(k1+1)ω0(t0+(N-1)Ts)
︙︙︙
e-jk2ω0t0e-jk2ω0(t0+Ts)…e-jk2ω0(t0+(N-1)Ts)x(t0)
x(t0+Ts)
︙
x(t0+(N-1)Ts)


（5.5）


矩阵向量乘法实现例5.1.1数值解的MATLAB源代码如下

close all; clc;clear all;

T = 4;%周期 

N = 400;%一个周期抽样点数

dt = T/N;%抽样间隔

t=-T/2:dt:T/2-dt;%时域自变量

x=2*rectpuls(t,2);%x(t)的赋值

X = ［］;%傅里叶级数

w0 = 2*pi/T;%基频

k = -10:10;%需要计算的Xk的项数

［W,tt］ = meshgrid(k*w0,t); %将自变量转为矩阵形式

Xk = dt/T*x*exp(-j*tt.*W);%利用矩阵-向量乘法计算Xk

stem(k*w0,Xk,'filled'); title('Xk')

程序运行结果如图5.1.4所示。在参数相同的前提下，采用for循环结构计算Xk的程序运行时间为0.010007s，而矩阵向量相乘计算Xk的程序运行时间为0.004057s。所以在MATLAB编程时需要尽量用矩阵形式实现算法。


图5.1.4矩阵向量乘法得到的连续时间周期矩形脉冲信号的频谱


观察图5.1.4可以发现，对于周期矩形脉冲信号的双边谱来说，幅度谱是偶函数，相位谱是奇函数。但不是所有周期信号的双边谱都有此特点，例如图5.1.5所示的虚指数信号ejπt2的双边谱。可以证明，只有实信号，才具有幅度谱是偶函数、相位谱是奇函数的特点。


图5.1.5ejπt2的双边谱


MATLAB提供了快速傅里叶变换(Fast Fourier Transform，FFT)函数fft来对离散时间信号进行频谱分析。可以先将连续时间信号抽样，再对得到的离散时间序列用fft函数计算傅里叶级数，近似得到连续时间信号的频谱，详见6.1.2节。
5.1.2周期信号的频谱
图5.1.2以角频率为横轴，Xk为纵轴画出了周期信号的双边谱，这是Xk为实数的情况。若Xk为复数，需要分别以Xk的模、初相为纵轴，得到的图形称为信号的幅度谱和相位谱，合称频谱图，借助频谱图可以实现对信号的频谱分析。

例5.1.2对如图5.1.1所示的周期矩形脉冲信号，分别改变信号的周期T和脉冲宽度τ，讨论其对频谱的影响。
解： 可以利用5.1.1节的方法，先得到连续时间周期信号在一个周期的抽样序列
x(n)，再通过数值计算得到
Xk。但由于周期矩形脉冲信号傅里叶级数Xk的理论值为


Xk=2τTSakω0τ2
（5.6）


故本题中直接根据理论公式计算Xk，然后对频谱进行分析。首先保持脉冲宽度不变，增加周期。MATLAB源代码如下

clc; clear all; close all;

dt = 1e-2; %抽样间隔

t=-20:dt:20;%时域自变量,画出3个完整周期

%保持脉冲宽度不变，改变周期

tao = 2;%脉冲宽度

T = 4:2:10;%周期

for i=1:4

x = square((t+tao/2)*(2*pi)/T(i),tao/T(i)*100)+1; %周期信号的赋值

w0 = 2*pi/T(i);%基频

k = -round(10/w0):round(10/w0);%画-10:10范围内的频谱

Xk = 2*tao/T(i)*sinc(k*w0*tao/2/pi);%Xk的理论值

subplot(4,2,2*i-1); plot(t,x); 

xlabel('t'); title('x(t)'); axis(［-20,20,-0.1,2.1］)

subplot(4,2,2*i);stem(k*w0,Xk,'filled');

xlabel('＼omega'); title('Xk'); axis(［-10 10 -0.2 1］)

end

程序运行结果如图5.1.6所示，动态结果请扫描二维码。















图5.1.6周期对频谱的影响



接下来，保持周期不变，增加脉冲宽度。MATLAB源代码如下

close all; clc;clear all;

dt = 1e-2;%抽样间隔

t=-20:dt:20;%时域自变量,画3个完整周期%保持周期不变，改变脉冲宽度

T = 15;%周期

for i = 1:4

tao = 2*i-1;%脉冲宽度

x = square((t+tao/2)*(2*pi)/T,tao/T*100)+1;%周期信号的赋值

w0 = 2*pi/T;%基频

k = -round(10/w0):round(10/w0);%画-10:10范围内的频谱

Xk = 2*tao/T*sinc(k*w0*tao/2/pi);%Xk的理论值

subplot(4,2,2*i-1); plot(t,x);

xlabel('t'); title('x(t)'); axis(［-20,20,-0.1,2.1］);

subplot(4,2,2*i); stem(k*w0,Xk,'filled');

xlabel('＼omega'); title('Xk'); axis(［-10,10,-0.2,1］)

hold on

plot(［-2*pi/tao,2*pi/tao］,［0,0］,'r','linewidth',2)%画主瓣宽度

hold off

end 

程序运行结果如图5.1.7所示，为便于观察频谱的主瓣宽度，图中用红色实线画出了频谱图中-2πτ,2πτ的范围。动态结果请扫描二维码。


























图5.1.7脉冲宽度对频谱的影响



观察发现： 
（1） 脉冲宽度τ保持不变，若T增大，则频谱主瓣高度EτT减小，各条谱线高度也相应地减小，各谱线间隔ω0=2πT减小，谱线变密； 若T减小，则情况相反。但不管T如何改变，频谱包络的第一个零点±2πτ不变，频谱主瓣宽度不变。
（2） 周期T保持不变，若τ增大，则频谱主瓣高度EτT增大，各条谱线高度也相应地增大，谱包络的第一个零点±2πτ减小，主瓣宽度变窄； 若τ减小，则情况相反。但不管τ如何改变，谱线间隔ω0=2πT不变。
综上所述，谱线间隔只与周期有关，且与其成反比； 频谱主瓣宽度仅与脉冲宽度有关，且与其成反比。而频谱高度与周期和脉冲宽度都有关，且与周期成反比，与脉冲宽度成正比。可以想象，当脉冲宽度保持不变，周期趋于无穷大时： 时域上，周期信号变成了非周期信号； 频域上，离散谱变成了连续谱，且高度为零，显然傅里叶级数不能用于非周期信号的频谱分析。
5.1.3周期信号的分解
周期信号三角形式傅里叶级数如下


x(t)=c0+∑∞k=1ckcos(kω0t+φk)
（5.7）


其中，ckcos(kω0t+φk)称为周期信号的第k次谐波，与式（5.2）相比，


c0=X0,ck=2|Xk|,φk=∠Xk
（5.8）


由周期信号得到各次谐波的过程称为信号的分解； 反过来，由各次谐波相加得到周期信号的过程称为信号的合成。如果需要直接观察真实存在的各次谐波，可以将周期信号通过滤波器组，详见11.1节周期矩形脉冲信号的分解与合成实验。但如果只是想观察各次谐波的特点，可以根据各次谐波理论上的振幅、初相，通过ckcos(kω0t+φk)仿真产生第k次谐波。
例5.1.3对如图5.1.1所示的周期矩形脉冲信号，画出其0~8次谐波。
解： 与例5.1.2类似，直接利用理论值计算得出Xk，再根据式（5.8）计算ck和φk，从而得出第k次谐波。MATLAB源代码如下

close all; clc;clear all;

tao = 2;%脉冲宽度

T = 4;%周期

dt = 1e-2;%抽样间隔

t=-T/2:dt:T/2;%时域自变量,取1个完整周期

w0 = 2*pi/T;%基频

k = 0:8;%谐波次数

Xk = 2*tao/T*sinc(k*w0*tao/2/pi);%Xk的理论值

Ck =［Xk(1),2*abs(Xk(2:end))］;%三角形式傅里叶级数振幅

faik = angle(Xk);

for i=1:9

yk=Ck(i)*cos(k(i)*w0*t+faik(i));

subplot(3,3,i); plot(t,yk);

xlabel('t'); title(［'第',num2str(k(i)),'次谐波'］);

axis(［-T/2,T/2,-1.5,1.5］)

end

程序运行结果如图5.1.8所示。


图5.1.8周期矩形脉冲信号的各次谐波


观察发现，偶次谐波均为0，这是因为该周期矩形脉冲信号减去其直流后为方波信号，属于半波像对称信号，半波像对称信号的偶次谐波理论上均为0； 奇次谐波角频率依次增加，振幅依次递减，符合周期信号频谱收敛性的特点。
5.1.4周期信号的合成
从式(5.7)可以看出，对于信号合成来说，需要无限次谐波相加才能逼近周期信号x(t)，但实际中只能实现有限次谐波相加。为观察不同谐波对信号合成的影响，用MATLAB仿真结果做进一步的说明。
例5.1.4对例5.1.3得到的各次谐波，分析各次谐波对信号合成的影响。
解： 依次增加参与合成的谐波次数，MATLAB源代码如下

close all; clc;clear all;

tao = 2;%脉冲宽度

T = 4;%周期

dt = 1e-2;%抽样间隔

t=-T/2:dt:T/2;%时域自变量,取1个完整周期

w0 = 2*pi/T;%基频

k = 0:8;%谐波次数

Xk = 2*tao/T*sinc(k*w0*tao/2/pi);%Xk的理论值

Ck =［Xk(1),2*abs(Xk(2:end))］;%三角形式傅里叶级数振幅

faik = angle(Xk);

x =2*rectpuls(t,tao);

for i=1:9

yk(i,:)=Ck(i)*cos(k(i)*w0*t+faik(i));

end

subplot(2,2,1); plot(t,x);

xlabel('t'); title('原信号'); ylim(［-0.5,2.5］);

subplot(2,2,2); plot(t,sum(yk(1:4,:),1));

xlabel('t'); title('前3次谐波参与合成'); ylim(［-0.5,2.5］);

subplot(2,2,3); plot(t,sum(yk(1:6,:),1));

xlabel('t'); title('前5次谐波参与合成'); ylim(［-0.5,2.5］);

subplot(2,2,4); plot(t,sum(yk(1:8,:),1));

xlabel('t'); title('前7次谐波参与合成'); ylim(［-0.5,2.5］);

程序运行结果如图5.1.9所示。观察发现，随着参与合成的谐波次数增加，合成波形越来越逼近周期信号，但间断点处总有一个过冲。可以证明，只要不是所有谐波参与合成，合成波形总会在间断点处有个过冲，且过冲值不变，这种现象即为Gibbs现象。


图5.1.9周期矩形脉冲信号合成的结果


为分析单次谐波对信号合成的影响，将前7次谐波的合成结果依次去除1、3、5、7次谐波，结果如图5.1.10所示。为便于比较，用点线画出了被分解的周期矩形脉冲信号。比较发现，低次谐波对信号合成的影响较大。这是因为低次谐波频率较小，反映了信号中缓慢变化的物理量； 高次谐波频率较大，反映了信号中快速变化的物理量。一般来说，信号中缓慢变化的部分占主导地位。



图5.1.10单次谐波对信号合成的影响


例5.1.5已知某周期三角波信号的周期为
4，其指数形式傅里叶系数为
Sa2(0.5kπ)，画出其前N次谐波合成的信号。
解： 不断增加参与合成的谐波次数，MATLAB源代码如下

close all; clc;clear all;

tao = 2;%脉冲宽度

T = 4;%周期

dt = 1e-2;%抽样间隔

t=-T/2:dt:T/2;%时域自变量,取1个完整周期

w0 = 2*pi/T;%基频

k = 0:11;%谐波次数

Xk = sinc(k/2).^2;%Xk的理论值

Ck =［Xk(1),2*abs(Xk(2:end))］;%三角形式傅里叶级数振幅

faik = angle(Xk);

for i=1:12

yk(i,:)=Ck(i)*cos(k(i)*w0*t+faik(i));

end

subplot(2,2,1); plot(t,sum(yk(1:4,:),1));

xlabel('t');title('前3次谐波参与合成')

subplot(2,2,2); plot(t,sum(yk(1:6,:),1));

xlabel('t');title('前5次谐波参与合成')

subplot(2,2,3); plot(t,sum(yk(1:10,:),1));

xlabel('t');title('前9次谐波参与合成')

subplot(2,2,4); plot(t,sum(yk(1:12,:),1));

xlabel('t');title('前11次谐波参与合成')

程序运行结果如图5.1.11所示。


图5.1.11周期三角波信号合成的结果


比较图5.1.9和图5.1.11可以发现，同等条件下显然周期三角波信号合成效果较好。这是因为三角波信号的各次谐波的复振幅是Sa2函数，而周期矩形脉冲信号的是Sa函数，而Sa2函数比Sa函数衰减快得多。正因为周期三角波信号各次谐波衰减更快，所以高次谐波对信号合成的影响相对更小。
5.2非周期信号的傅里叶变换
为有效分析非周期信号的频谱，引入了傅里叶变换。傅里叶变换定义为


X(jω)=∫∞-∞x(t)e-jωtdt
（5.9）


傅里叶反变换定义为


x(t)=12π∫∞-∞X(jω)ejωtdω
（5.10）


x(t)与X(jω)一一对应，构成傅里叶变换对，记作

x(t)FX(jω)

5.2.1傅里叶变换的符号函数求解
MATLAB符号工具箱提供了求解傅里叶变换的函数fourier和求解傅里叶反变换的函数ifourier，它们的调用格式如下。
fourier(x)： 对默认变量为t的符号表达式求傅里叶变换，默认返回关于w的函数。
fourier(x,v)： 对默认自变量为t的符号表达式求傅里叶变换，返回关于v的函数。
fourier(x,u,v)： 对x(u)求傅里叶变换，返回关于v的函数。
ifourier(X)： 对默认变量为w的符号函数表达式求傅里叶反变换，默认返回关于t的函数。
ifourier(X,u)： 对默认变量为w的符号函数表达式求傅里叶反变换，返回关于u的函数。
ifourier(X,v,u)： 对X(v)求傅里叶反变换，返回关于u的函数。
需要注意的是，在调用函数fourier和ifourier之前，要用syms函数定义所用到的符号变量或符号表达式，返回结果中默认用i表示虚数单位。
可利用fplot实现对符号函数的画图，调用格式如下。
fplot(f)： 在默认区间［-5 5］，绘制由函数f定义的曲线。
fplot(f,xinterval)： 在xinterval指定的区间绘图，xinterval是［xmin xmax］形式的二元素向量。
例5.2.1求e-tu(t)的傅里叶变换，并画出其频谱图。

解： MATLAB源代码如下

close all; clc;clear all;

syms t w

Xjw = fourier(exp(-t)*heaviside(t),w) %用符号函数求傅里叶变换

subplot(1,2,1); fplot(abs(Xjw),［-10 10］); title('e^{-t}u(t)的幅度谱')

subplot(1,2,2); fplot(angle(Xjw),［-10 10］); title('e^{-t}u(t)的相位谱')

程序运行结果如图5.2.1所示。


图5.2.1e-tu(t)的傅里叶变换


命令窗口运行结果为

Xjw = 

1/(1 + w*1i)


例5.2.2求1ω2+1的傅里叶反变换，并画出波形图。
解： MATLAB源代码如下

close all; clc;clear all;

syms t w

xt = ifourier(1/(w^2+1),w,t)%用符号函数求傅里叶变换

fplot(xt); ylim(［0 0.5］); title('1/(＼omega^{2}+1)的傅里叶反变换')

程序运行结果如图5.2.2所示。


图5.2.21ω2+1的傅里叶反变换


命令窗口运行结果为

xt = 

exp(-abs(t))/2



5.2.2傅里叶变换的数值计算
当信号的数学表达式非常复杂，或者得到的是信号的抽样结果、无法写出其数学表达式时，此时无法用fourier函数得出信号的傅里叶变换，但可以借助数值计算的方法进行求解。 

下面讨论利用数值求解计算连续时间非周期信号傅里叶变换的方法。假设在非周期信号的主要取值区间［t1,t2］内抽样了N个点，则抽样间隔Ts=t2-t1N，与式（5.3）类似，有


X(jω)≈Ts∑N-1n=0x(t1+nTs)e-jω(t1+nTs)
（5.11）


用式（5.11）可以计算出任意频点的傅里叶变换值。假设非周期信号频谱的主要取值区间为［ω1,ω2］，在其间均匀抽样了M个值，则



Δω=ω2-ω1M 

X［j(ω1+mΔω)］≈Ts∑N-1n=0x(t1+nTs)e-j(ω1+mΔω)(t1+nTs)
（5.12）


可以采用同样的方法计算傅里叶反变换


x(t1+nTs)≈Δω2π∑M-1m=0X［j(ω1+mΔω)］ej(ω1+mΔω)(t1+nTs)
（5.13）


式（5.12）、式（5.13）可用矩阵向量乘法编程实现，以式（5.12）为例


X(jω1)
X［j(ω1+Δω)］
︙
X［j(ω2-Δω)］≈Tse-jω1t1e-jω1(t1+Ts)…e-jω1(t2-Ts)
e-j(ω1+Δω)t1e-j(ω1+Δω)(t1+Ts)…e-j(ω1+Δω)(t2-Ts)
︙︙︙
e-j(ω2-Δω)t1e-j(ω2-Δω)(t1+Ts)…e-j(ω2-Δω)(t2-Ts)x(t1)
x(t1+Ts)
︙
x(t2-Ts)


（5.14）


例5.2.3用数值计算的方法求e-tu(t)的傅里叶变换。
解： 采用矩阵向量相乘实现傅里叶变换的数值解，MATLAB源代码如下

close all; clc;clear all;

dt = 0.01;%时域抽样间隔

t=0:dt:20;%信号主要取值区间 

x=exp(-t); %信号赋值

w = -20:0.01:20; %信号频谱主要取值区间

［W,T］ = meshgrid(w,t); %生成矩阵

X = dt*x*exp(-j*T.*W); %利用矩阵-向量乘法计算

subplot(1,2,1); plot(w,abs(X)); title('e^{-t}u(t)的幅度谱')

subplot(1,2,2); plot(w,angle(X)); title('e^{-t}u(t)的相位谱')

程序运行结果如图5.2.3所示。


图5.2.3矩阵向量乘法得到的傅里叶变换


将矩阵向量乘法得到的数值解与理论值1jω+1进行比较，结果如图5.2.4所示。可以看出计算误差略大，增加t的范围，结果几乎没有改善，这是因为前面已经取到了信号的主要取值区间； 将时域抽样间隔由10-2降低到10-4，误差将降到10-5量级，但运算量会显著增加。


图5.2.4矩阵向量乘法的误差


在抽样间隔相同的条件下，5.1.1节中矩阵向量乘法计算傅里叶系数的误差却是10-4量级，一方面是因为积分区间只需取一个完整的周期，不需要取(-∞，∞)； 另一方面是因为Xk是离散的，只需计算有限个频率点上的复振幅。

MATLAB提供了quad8和quadl函数来计算数值积分，quad8函数的返回值是用自适应Simpson算法得出的积分值。quadl函数是从MATLAB 6.0版本才开始出现的一个积分函数，它的返回值是用Lobatto算法得到的积分值，具有更高的积分精度。quadl的调用格式为

y = quadl(fun,a,b)

其中,a,b分别是积分下限和积分上限，fun是被积函数。
用quadl函数重新求解例5.2.3，MATLAB源代码如下

close all; clc;clear all;

w = -20:0.01:20; %信号频谱主要取值区间

for i=1:length(w)

F = @(t)exp(-t).*exp(-j*w(i)*t); %被积函数

X(i)=quadl(F,0,20);%Lobatto法

end

figure;plot(abs(X)-abs(1./(j*w+1)));title('幅度谱计算误差')

程序运行结果如图5.2.5所示。比较发现，Lobatto算法误差比矩阵向量乘法的误差小得多。但这种方法需要知道信号的函数表达式，如果只能得到信号的抽样结果，仍旧需要用矩阵向量乘法计算傅里叶变换。


图5.2.5Lobatto算法计算数值积分的误差

5.3傅里叶变换性质
傅里叶变换性质是信号与系统课程的重要知识点，利用傅里叶变换的性质不仅可以简化计算，而且可以更好地理解时域特性与频域特性间的关系以及时域变化与频域变化间的关系。
5.3.1奇偶特性
理论上，实偶信号的频谱为实偶函数，实奇信号的频谱为虚奇函数。但这并不意味着实偶信号的幅度谱、相位谱均是偶函数，实奇信号的幅度谱、相位谱均是奇函数。频谱是实偶函数指的是频谱是实的，且


X(jω)=X(-jω)
（5.15）


将X(jω)写成直角坐标形式，可以推出实偶信号频谱的实部、虚部均是偶函数，


Re［X(jω)］=Re［X(-jω)］
（5.16）

Im［X(jω)］=Im［X(-jω)］


将X(jω)写成极坐标形式，可以推出幅度谱是偶函数，且ej∠X(jω)=ej∠X(-jω)，但并不能由此得出相位谱是偶函数的结论。
同理，实奇信号的频谱是虚奇函数，指的是频谱是虚的，且频谱的实部、虚部都是奇函数。下面通过例5.3.1验证以上分析。
例5.3.1分别画出G4(t)、Λ4(t)、e-tu(t)、e-tu(t)-etu(-t)的频谱。
解： 分别画出4个信号的时域波形、幅度谱和相位谱，MATLAB源代码如下

close all; clc;clear all;

dt = 1e-3;%时域抽样间隔

t=-10:dt:10;%信号主要取值区间 

w = -20:0.01:20; %信号频谱主要取值区间

［W,T］ = meshgrid(w,t); %生成矩阵

%门信号

tao=4;

xt1 = rectpuls(t,tao);

Xjw1 = dt*xt1*exp(-j*T.*W); %利用矩阵-向量乘法计算

%Sa(t)

xt2 = tripuls(t,tao);

Xjw2 = dt*xt2*exp(-j*T.*W); %利用矩阵-向量乘法计算

%e^(-t)u(t)

xt3 = exp(-t).*(t>=0);

Xjw3 = dt*xt3*exp(-j*T.*W); %利用矩阵-向量乘法计算

%sgn(t)

xt4 = exp(-t).*(t>=0)-exp(t).*(t<=0);

Xjw4 = dt*xt4*exp(-j*T.*W); %利用矩阵-向量乘法计算

%画图

subplot(4,3,1); plot(t,xt1);

title('G_{4}(t)'); xlabel('t'); ylim(［-0.1 1.1］)

subplot(4,3,2); plot(w,abs(Xjw1));

title('G_{4}(t)的幅度谱'); xlabel('＼omega')

subplot(4,3,3); plot(w,angle(Xjw1).*(abs(Xjw1)>=1e-3));%去除数值计算带来的误差

title('G_{4}(t)的相位谱'); xlabel('＼omega'); ylim(［-3.2 3.2］)

subplot(4,3,4); plot(t,xt2);

title('＼Lambda_{4}(t)'); xlabel('t')

subplot(4,3,5); plot(w,abs(Xjw2));

title('＼Lambda_{4}(t)的幅度谱'); xlabel('＼omega')

subplot(4,3,6); plot(w,angle(Xjw2).*(abs(Xjw2)>=1e-3));%去除数值计算带来的误差

title('＼Lambda_{4}(t)的相位谱'); xlabel('＼omega'); ylim(［-3.2 3.2］)

subplot(4,3,7); plot(t,xt3);

title('e^{-t}u(t)'); xlabel('t'); ylim(［-0.1 1.1］)

subplot(4,3,8); plot(w,abs(Xjw3));

title('e^{-t}u(t)的幅度谱'); xlabel('＼omega')

subplot(4,3,9); plot(w,angle(Xjw3).*(abs(Xjw3)>=1e-3));%去除数值计算带来的误差

title('e^{-t}u(t)的相位谱'); xlabel('＼omega'); ylim(［-3.2 3.2］)

subplot(4,3,10); plot(t,xt4);

title('e^{-t}u(t)-e^{t}u(-t)'); xlabel('t'); ylim(［-1.1 1.1］)

subplot(4,3,11); plot(w,abs(Xjw4));

title('e^{-t}u(t)-e^{t}u(-t)的幅度谱'); xlabel('＼omega')

subplot(4,3,12); plot(w,angle(Xjw4).*(abs(Xjw4)>=1e-3));%去除数值计算带来的误差

title('e^{-t}u(t)-e^{t}u(-t)的相位谱'); xlabel('＼omega'); ylim(［-3.2 3.2］) 

程序运行结果如图5.3.1所示。在画相位谱时，考虑到幅度较小时容易带来相位计算误差，因此认为当幅度谱小于10-3时，相位为0。从图5.3.1中可以看出，虽然G4(t)、Λ4(t)均为实偶信号，但相位谱并不均是偶函数，与前文的分析一致； 同样的，虽然e-tu(t)-etu(-t)是实奇信号，但其幅度谱也不是奇函数。4个信号都是实信号，满足实信号幅度谱均是偶函数，相位谱均是奇函数的特点。 


图5.3.14个信号的时频图


分别画出4个信号频谱的实部和虚部，结果如图5.3.2所示。与前面的分析一致，偶信号频谱的实部是偶函数，虚部为零； 奇信号频谱的虚部是奇函数，实部为零； 既非奇又非偶的信号频谱实部、虚部均不全为零，但也满足实信号频谱的实部是偶函数、虚部是奇函数的特点。


图5.3.24个信号频谱的实部、虚部


5.3.2展缩特性
若x(t)FX(jω)，由展缩性质可知


x(at)F1|a|Xjωa
（5.17）


这意味着信号在时域上压缩，频域上将扩展； 反过来在时域上扩展，频域上将压缩。下面通过例5.3.2说明这一点。
例5.3.2设x(t)=G4(t)，分别画出x(t)、x(t/2)、x(2t)的频谱。
解： MATLAB源代码如下


close all; clc;clear all;

dt = 1e-3;%时域抽样间隔

t=-10:dt:10;%信号主要取值区间 

w = -20:0.01:20; %信号频谱主要取值区间

［W,T］ = meshgrid(w,t); %生成矩阵

%x(t)

tao=4;

xt1 = rectpuls(t,tao);

Xjw1 = dt*xt1*exp(-j*T.*W); %利用矩阵-向量乘法计算

%x(t/2)

xt2 = rectpuls(t/2,tao);

Xjw2 = dt*xt2*exp(-j*T.*W); %利用矩阵-向量乘法计算

%x(2t)

xt3 = rectpuls(2*t,tao);

Xjw3 = dt*xt3*exp(-j*T.*W); %利用矩阵-向量乘法计算

%画图

subplot(3,3,1); plot(t,xt1);

title('x(t)'); xlabel('t'); ylim(［-0.1 1.1］)

subplot(3,3,2); plot(w,abs(Xjw1));

title('x(t)的幅度谱'); xlabel('＼omega')

subplot(3,3,3); plot(w,angle(Xjw1).*(abs(Xjw1)>=1e-3));%去除数值计算带来的误差

title('x(t)的相位谱'); xlabel('＼omega');ylim(［-3.2 3.2］)

subplot(3,3,4); plot(t,xt2);

title('x(t/2)'); xlabel('t'); ylim(［-0.1 1.1］)

subplot(3,3,5); plot(w,abs(Xjw2));

title('x(t/2)的幅度谱'); xlabel('＼omega')

subplot(3,3,6); plot(w,angle(Xjw2).*(abs(Xjw2)>=1e-3));%去除数值计算带来的误差

title('x(t/2)的相位谱'); xlabel('＼omega'); ylim(［-3.2 3.2］)

subplot(3,3,7); plot(t,xt3);

title('x(2t)'); xlabel('t'); ylim(［-0.1 1.1］)

subplot(3,3,8); plot(w,abs(Xjw3));

title('x(2t)的幅度谱'); xlabel('＼omega')

subplot(3,3,9); plot(w,angle(Xjw3).*(abs(Xjw3)>=1e-3));%去除数值计算带来的误差

title('x(2t)的相位谱'); xlabel('＼omega'); ylim(［-3.2 3.2］)


程序运行结果如图5.3.3所示。可以看出，从x(t)到x(t/2)，时域上的宽度扩展为原来的2倍，频带宽度却压缩为原来的一半，同时幅度谱最大值变成原来的2倍，从理论上说，式（5.17）中有一个1|a|的系数，也可以直观理解为信号在时域上变宽，本身的能量变大，而幅度谱却变窄，结果只能是幅度谱的值变大。另一方面，从x(t)到x(2t)，时域上的宽度压缩为原来的一半，频带宽度却扩展为原来的2倍，同时幅度谱最大值变成原来的一半。这些变化直观反映了傅里叶变换的展缩特性。


图5.3.3信号展缩前后的频谱







可以通过语音信号的展缩更加直观地感受时域展缩对语速、音调的影响，请扫描二维码。
5.3.3时移特性
若x(t)FX(jω)，由时移性质可知


x(t+t0)FX(jω)ejωt0
（5.18）


这意味着信号时移后，频谱将乘以因子ejωt0。由于两复数相乘对应模相乘、相位相加，再加上ejωt0的模为1，相位为ωt0，因此，信号的时移不会影响幅度谱，只是相位谱产生附加变化ωt0，改变量是过原点且斜率是t0的直线。下面通过例5.3.3说明这一点。
例5.3.3分别画出x(t)=Λ4(t)、x(t-0.1)、x(t-1)的频谱。
解： 为方便比较，选用相位谱为0的Λ4(t)作为x(t)。MATLAB源代码如下

close all; clc;clear all;

dt = 1e-2;%时域抽样间隔

t=-10:dt:10;%信号主要取值区间 

w = -10:0.01:10; %信号频谱主要取值区间

［W,T］ = meshgrid(w,t); %生成矩阵

%x(t)

tao=4;

xt1 = tripuls(t,tao);

Xjw1 = dt*xt1*exp(-j*T.*W); %利用矩阵-向量乘法计算

%x(t-0.1)

xt2 = tripuls(t-0.1,tao);

Xjw2 = dt*xt2*exp(-j*T.*W); %利用矩阵-向量乘法计算

%x(t-1)

xt3 = tripuls(t-1,tao);

Xjw3 = dt*xt3*exp(-j*T.*W); %利用矩阵-向量乘法计算

%画图

subplot(3,3,1); plot(t,xt1);

title('x(t)'); xlabel('t'); ylim(［-0.1 1.1］)

subplot(3,3,2); plot(w,abs(Xjw1));

title('x(t)的幅度谱'); xlabel('＼omega')

subplot(3,3,3); plot(w,angle(Xjw1));

title('x(t)的相位谱'); xlabel('＼omega');ylim(［-3.2 3.2］)

subplot(3,3,4); plot(t,xt2);

title('x(t-0.1)'); xlabel('t'); ylim(［-0.1 1.1］)

subplot(3,3,5); plot(w,abs(Xjw2));

title('x(t-0.1)的幅度谱'); xlabel('＼omega')

subplot(3,3,6); plot(w,angle(Xjw2));

title('x(t-0.1)的相位谱'); xlabel('＼omega'); ylim(［-3.2 3.2］)

subplot(3,3,7); plot(t,xt3);

title('x(t-1)'); xlabel('t'); ylim(［-0.1 1.1］)

subplot(3,3,8); plot(w,abs(Xjw3));

title('x(t-1)的幅度谱'); xlabel('＼omega')

subplot(3,3,9); plot(w,angle(Xjw3));

title('x(t-1)的相位谱'); xlabel('＼omega'); ylim(［-3.2 3.2］)


程序运行结果如图5.3.4所示。可以看出，信号时移前后幅度谱没有发生改变，只是相位谱发生了变化，相位改变量是过原点的直线，该直线的斜率跟时移量有关，结合式（5.18），直线的斜率刚好就是时移量。x(t-1)的相位谱之所以是分段直线而非前面分析的直线，是因为要保证初相在［-π,π］区间内。


图5.3.4信号时移前后的频谱


5.3.4频移特性
若x(t)FX(jω)，由频移性质可知


x(t)ejω0tFX［j(ω-ω0)］
（5.19）


这意味着信号乘以虚指数信号ejω0t后，频谱将移位ω0。若信号乘以余弦信号cosω0t，通过将cosω0t进行欧拉公式展开，可得


x(t)cosω0tFX［j(ω-ω0)］+X［j(ω+ω0)］2
（5.20）


这意味着信号乘以cosω0t后，频谱将保持形状不变，往左往右各频移ω0，同时幅度将变成原来的一半。下面通过例5.3.4说明这一点。
例5.3.4分别画出G4(t)、G4(t)cos(50t)的频谱。
解： MATLAB源代码如下

close all; clc;clear all;

dt = 1e-2;%时域抽样间隔

t=-10:dt:10;%信号主要取值区间 

w = -60:0.01:60; %信号频谱主要取值区间

［W,T］ = meshgrid(w,t); %生成矩阵

%x(t)

tao=4;

xt1 = rectpuls(t,tao);

Xjw1 = dt*xt1*exp(-j*T.*W); %利用矩阵-向量乘法计算

%x(t)cos(50t)

xt2 = rectpuls(t,tao).*cos(50*t);

Xjw2 = dt*xt2*exp(-j*T.*W); %利用矩阵-向量乘法计算

%画图

subplot(2,2,1); plot(t,xt1);

title('G_{4}(t)'); xlabel('t'); ylim(［-0.1 1.1］)

subplot(2,2,2); plot(w,Xjw1);

title('G_{4}(t)的频谱'); xlabel('＼omega')

subplot(2,2,3); plot(t,xt2);

title('G_{4}(t)cos(50t)'); xlabel('t')

subplot(2,2,4); plot(w,Xjw2);

title('G_{4}(t)cos(50t)的频谱'); xlabel('＼omega')

程序运行结果如图5.3.5所示。可以看出，乘以cos(50t)后，G4(t)频谱的中心点由0移到了±50，同时最大值由4降为2，与前面的理论分析一致。


图5.3.5频移前后信号的频谱



5.4连续时间系统的频率响应
对式(4.17)描述的连续时间线性时不变系统等式两边分别做傅里叶变换并进行整理，可得


H(jω)=bm(jω)m+bm-1(jω)m-1+…+b0(jω)n+an-1(jω)n-1+…+a0
（5.21）


式中，m和n都是正整数，系数均为实数。
在MATLAB中，信号处理工具箱中的freqs函数可直接计算连续时间系统的频率响应，调用格式如下。
freqs(b,a)： 没有返回值，直接画出系统幅频、相频响应的波特图。b、a分别是式(5.21)中的分子、分母多项式的系数向量。
H = freqs(b,a,w)： 向量w是系统频率响应的角频率范围，返回值H为w上的频率响应。
［H,w］ = freqs(b,a,n)： 自动设定n个角频率点来计算频率响应，返回值w为设定的n个角频率值，n的默认值为200。
例5.4.1某系统微分方程为y″(t)+y′(t)+y(t)=x(t)，求该系统的频率响应并画出幅频、相频响应曲线。
解： MATLAB源代码如下

close all; clc;clear all;

a = ［1 1 1］;%分母多项式系数

b = ［1］;%分子多项式系数

［H,w］=freqs(b,a);%求连续时间系统的频响

subplot(1,2,1); plot(w,abs(H)); 

xlabel('＼omega'); title('幅频响应')

subplot(1,2,2); plot(w,angle(H)); 

xlabel('＼omega'); title('相频响应')

程序运行结果如图5.4.1所示。


图5.4.1连续时间系统的频率响应


调用freqs函数时采用了自动设定自变量的方式，结果只画出了频谱的右半边。可以通过设置自变量，得到双边谱。将上述程序改为以下形式，得到的双边谱如图5.4.2所示。

close all; clc;clear all;

a = ［1 1 1］;%分母多项式系数

b = ［1］;%分子多项式系数

w = -10:0.01:10;%角频率

H =freqs(b,a,w);%求连续时间系统的频响

subplot(1,2,1); plot(w,abs(H));

xlabel('＼omega'); title('幅频响应')

subplot(1,2,2); plot(w,angle(H));

xlabel('＼omega'); title('相频响应')




图5.4.2程序改进后连续时间系统的频率响应


5.5正弦信号通过系统的响应
正余弦信号通过系统的响应既可以通过2.1.2节零状态响应的时域求解来实现，也可以快速地通过式(5.22)、式(5.23)以及叠加原理得到理论结果。


cos(ω0t)→|H(jω0)|cos(ω0t+∠H(jω0))
（5.22）

sin(ω0t)→|H(jω0)|sin(ω0t+∠H(jω0))
（5.23）


例5.5.1计算x(t)=cost+cos10t通过例5.4.1系统的响应。
解： MATLAB源代码如下

clc; clear all; close all;

a=［1 1 1］;b=［1］; %系统分母、分子多项式系数

H = freqs(b,a,［1 10］); %仅计算在1、10两个角频率上的频率响应

t = 0:0.001:40;

x = cos(t)+cos(10*t);

y = abs(H(1))*cos(t+angle(H(1)))+abs(H(2))*cos(10*t+angle(H(2)));

subplot(1,2,1); plot(t,x);

xlabel('t'); title('x(t)')

subplot(1,2,2); plot(t,y);

xlabel('t'); title('y(t)')

程序运行结果如图5.5.1所示。观察发现，信号通过系统后高频分量几乎消失，这是因为系统的幅频响应在ω=10处几乎为0，而在ω=1处的值较大（也可以观察图5.4.2）。


图5.5.1余弦信号通过系统的响应

5.6无失真传输
5.6.1无失真传输系统

若某连续时间系统对所有的输入信号都可以实现


yzs(t)=kx(t-td)
（5.24）


即输出信号的幅度是输入信号幅度的k倍(k为实常数，且k≠0)，输出信号比输入信号延迟了td（td≥0），那么称该系统为无失真传输系统。对式（5.24）进行推导，得出无失真传输系统的幅频特性、相频特性须满足


|H(jω)|=k
（5.25）

∠H(jω)=-ωtd
（5.26）


例5.6.1某系统的微分方程为y′(t)+y(t)=x′(t)-x(t)，判断该系统是否是无失真传输系统。
解： 直接画出该系统的频率响应，MATLAB源代码如下

clc; clear all; close all;

a = ［1 1］;%分母多项式系数

b = ［1 -1］;%分子多项式系数

freqs(b,a);%直接画出连续时间系统的频响

程序运行结果如图5.6.1所示。观察发现，虽然系统的幅频特性恒为1，满足式(5.25)，但是相频特性不是直线，不满足式(5.26)，因此该系统不是无失真传输系统。


图5.6.1例5.6.1系统的频率响应


5.6.2信号的无失真传输
对于信号来说，只要在其频带范围内，系统的频率响应满足式(5.25)、式(5.26)，该信号通过系统后是无失真的。
例5.6.2某系统频率响应如图5.6.2所示，分别判断信号Sa(5t)、Sa(5t)ej5t、G5(t)、cos(6t)通过该系统后是否失真。


图5.6.2例5.6.2系统的频率响应


解： 利用符号函数分别计算4个信号的傅里叶变换，MATLAB源代码如下

close all; clc;clear all;

syms t w;%声明符号变量

%用符号函数求傅里叶变换

Xjw1 = fourier(sinc(5/pi*t),w); 

Xjw2 = fourier(sinc(5/pi*t)*exp(j*5*t),w); 

Xjw3 = fourier(heaviside(t+2.5)-heaviside(t-2.5),w);

Xjw4 = fourier(cos(6*t),w)

subplot(2,2,1); fplot(abs(Xjw1),［-12 12］); 

xlabel('＼omega');ylim(［-0.1 0.7］); title('Sa(5t)的幅度谱')

subplot(2,2,2); fplot(abs(Xjw2),［-12 12］);

xlabel('＼omega'); ylim(［-0.1 0.7］); title('Sa(5t)e^{j5t}的幅度谱')

subplot(2,2,3); fplot(abs(Xjw3),［-12 12］);

xlabel('＼omega'); ylim(［-0.1 5.1］); title('G_{5}(t)的幅度谱')

subplot(2,2,4); fplot(abs(Xjw4),［-12 12］); 

xlabel('＼omega'); title('cos(6t)的幅度谱')

程序运行结果如图5.6.3所示。观察发现： 
（1） Sa(5t)频带范围是［-5 5］，理论上通过系统后不失真； 
（2） Sa(5t)ej5t频带范围是［0 10］，虽然系统幅频响应在这个范围内是常数2，但相频响应不是直线，理论上通过系统后将出现失真； 
（3） G5(t)的频带范围超过了［-5 5］，理论上通过系统后也将出现失真； 
（4） cos(6t)的幅度谱显示的是零，但命令窗口的输出是Xjw4=pi*(dirac(w-6)+dirac(w+6))，也就是说，cos(6t)的频谱在±6上是无穷大，只是无法在图中显示出来。


图5.6.34个信号的幅度谱


接下来从理论上分析cos(6t)通过系统后是否失真。非常容易出错的一点是有读者认为cos(6t)的频谱在±6处有非零值，故其频带范围为［-6，6］。但频带范围指的是频谱非0的频率点的集合，而cos(6t)只有-6、6两个角频率，在这两个角频率上，幅度谱是常数，相位谱是过原点且斜率是-56的直线，因此该信号通过系统后无失真。或者由
式（5.22）可知


cos(6t)→|H(j6)|cos(6t+∠H(j6))=2cos(6t-5)
（5.27）


满足式（5.24）的定义。

例5.6.3画出例5.6.2中各信号通过系统的响应，验证该例题的结论。
解： 按照先将输入信号的傅里叶变换与系统的频率响应相乘，再傅里叶反变换求输出的思路，MATLAB源代码如下

close all; clc;clear all;

dw = 1e-2; %频谱抽样间隔

w = -12:dw:12; %信号频谱主要取值区间

%表示H(jw)

Hjw_angle = -5*sign(w).*(1-rectpuls(w,10))-w.*rectpuls(w,10);

Hjw = 2*rectpuls(w,20).*exp(j*Hjw_angle);

%考虑到cos信号，用矩阵-向量乘法计算X(jw)

dt = 1e-2;%时域抽样间隔

t=-10:dt:10;%信号主要取值区间 

xt1=sinc(5/pi*t); %信号赋值

xt2 = sinc(5/pi*t).*exp(j*5*t);

xt3 = rectpuls(t,5);

xt4 = cos(6*t);

［W,T］ = meshgrid(w,t); %生成矩阵

WT = exp(-j*T.*W); %傅里叶变换的矩阵

Xjw1 = dt*xt1*WT; %利用矩阵-向量乘法计算

Xjw2 = dt*xt2*WT; %利用矩阵-向量乘法计算

Xjw3 = dt*xt3*WT; %利用矩阵-向量乘法计算

Xjw4 = dt*xt4*WT; %利用矩阵-向量乘法计算

%输出傅里叶变换

Yjw1 = Hjw.*Xjw1;

Yjw2 = Hjw.*Xjw2;

Yjw3 = Hjw.*Xjw3;

Yjw4 = Hjw.*Xjw4;

%矩阵-向量法求傅里叶反变换

yt1 = dw/(2*pi)*Yjw1*WT';

yt2 = dw/(2*pi)*Yjw2*WT';

yt3 = dw/(2*pi)*Yjw3*WT';

yt4 = dw/(2*pi)*Yjw4*WT';

subplot(4,2,1); plot(t,xt1);

title('Sa(5t)'); xlim(［-5 5］)；

subplot(4,2,2); plot(t,yt1);

title('Sa(5t)通过系统后的输出'); xlim(［-5 5］)

subplot(4,2,3); plot(t,real(xt2)); 

title('Sa(5t)e^{j5t}的实部'); xlim(［-5 5］)

subplot(4,2,4); plot(t,real(yt2));

title('Sa(5t)e^{j5t}通过系统后输出的实部'); xlim(［-5 5］)

subplot(4,2,5); plot(t,xt3);

title('G_{5}(t)'); xlim(［-5 5］); ylim(［-0.1 1.1］)

subplot(4,2,6); plot(t,yt3);

title('G_{5}(t)通过系统后的输出'); xlim(［-5 5］); ylim(［-0.5 2.8］)

subplot(4,2,7); plot(t,xt4);

title('cos(6t)'); xlim(［-5 5］)

subplot(4,2,8); plot(t,yt4);

title('cos(6t)通过系统后的输出'); xlim(［-5 5］)

程序运行结果如图5.6.4所示。观察发现Sa(5t)、cos(6t)通过系统后，输出信号波形除幅度等比例增大以及有一些延时外，未发生其他变化，信号
通过系统后
未出现失真； 而Sa(5t)ej5t、G5(t)的波形发生了改变，信号通过系统后出现了失真，例5.6.2的结论正确。


图5.6.4例5.6.3输入、输出信号波形

5.7理想与实际滤波器
5.7.1理想滤波器

例5.7.1理想低通滤波器的频率响应如图5.7.1所示。分别讨论截止角频率ωc和td对滤波器单位冲激响应的影响。


图5.7.1理想低通滤波器的频率响应


解： 由于可以写出频率响应的解析表达式，所以可以利用符号函数计算频率响应的傅里叶反变换，得到单位冲激响应。先保持相频响应的斜率-td不变，增加ωc。MATLAB源代码如下

close all; clc;clear all;

syms t w

td = 1;%相频特性的负斜率

wc = 1:3;

for i=1:3

Hjw = (heaviside(w+wc(i))-heaviside(w-wc(i)))*exp(-j*w*td);

xt = ifourier(Hjw,w,t);%用符号函数求傅里叶变换

subplot(3,3,3*i-2); fplot(abs(Hjw)); 

ylim(［-0.1 1.1］);xlabel('＼omega');

title(［'＼omega_c=',num2str(wc(i)),'时的幅频响应'］)

subplot(3,3,3*i-1); fplot(angle(Hjw)); xlabel('＼omega'); 

title(［'＼omega_c=',num2str(wc(i)),'时的相频响应'］)

subplot(3,3,3*i); fplot(xt); 

ylim(［-0.3 1.3］); xlabel('t');

title(［'＼omega_c=',num2str(wc(i)),'时的单位冲激响应'］)

xticks(［-5 0 td 5］)%在横坐标上显示td 

end

程序运行结果如图5.7.2所示，动态结果请扫描二维码。随着截止角频率ωc的增加，单位冲激响应保持中心点在td处不变，但脉冲宽度不断减小，最大值不断增加。







图5.7.2ωc对滤波器单位冲激响应的影响


再保持ωc不变，增加相频响应的负斜率td。MATLAB源代码如下

close all; clc;clear all;

syms t w

wc = 2;

td = 0.5:1:2.5;

for i=1:3

Hjw = (heaviside(w+wc)-heaviside(w-wc))*exp(-j*w*td(i));

xt = ifourier(Hjw,w,t);%用符号函数求傅里叶变换

subplot(3,3,3*i-2); fplot(abs(Hjw)); 

ylim(［-0.1 1.1］); xlabel('＼omega');

title(［'td=',num2str(td(i)),'时的幅频响应'］)

subplot(3,3,3*i-1); fplot(angle(Hjw));

ylim(［-pi pi］); xlabel('＼omega'); 

title(［'td=',num2str(td(i)),'时的相频响应'］)

subplot(3,3,3*i); fplot(xt); xlabel('t');

title(［'td=',num2str(td(i)),'时的单位冲激响应'］)

xticks(［-5 0 td(i) 5］)%在横坐标上显示td 

end

程序运行结果如图5.7.3所示，动态结果请扫描二维码。随着td的增加，单位冲激响应保持脉冲宽度、最大值不变，中心不断右移。单位冲激响应的中心点即为相频响应的负斜率td。







图5.7.3td对滤波器单位冲激响应的影响


5.7.2实际滤波器
若根据系统的微分方程来判断该系统的滤波特性，可通过画出系统的频率响应实现。
例5.7.2某因果系统的微分方程为y″(t)+y′(t)+y(t)=x″(t)，判断该系统的滤波特性。
解： MATLAB源代码如下

close all; clc;clear all;

a = ［1 1 1］;%分母多项式系数

b = ［1 0 0］;%分子多项式系数

freqs(b,a);%画出频谱图的波特图

程序运行结果如图5.7.4所示，观察幅度特性发现该系统为高通滤波器。


图5.7.4某系统的频率响应







5.7.3利用Simulink实现信号的滤波
Simulink中提供了滤波器模块用以实现信号的滤波。这些模块位于DSP System Toolbox＼Filtering＼Filter Designs子库中。
例5.7.3借助Simulink，将信号x(t)=sin(2π·500t)+sin(2π·1000t)+sin(2π·2500t)中的每个正弦信号分别取出来。

解： Simulink模型如图5.7.5所示。通过将3个正弦信号相加产生需要的输入信号x(t)。需注意的是，在用Sine Wave模块产生正弦信号时，Sample time要小于1/(2×2500)，即抽样频率要大于5000Hz，并且该频率要与后面的滤波器模块中的抽样频率保持一致，才可以得到正确的滤波结果。这里抽样频率设置为44100Hz，该频率是常用的声音抽样频率。


图5.7.5例5.7.3的Simulink模型图


分别采用Lowpass、Bandpass、Highpass模块进行滤波，以得到sin(2π·500t)、sin(2π·1000t)、sin(2π·2500t)三个信号，注意滤波器模块中频率参数要留有一定的余量。根据需求分析，低通滤波器要保证通带截止频率大于500Hz，阻带截止频率小于1000Hz； 高通滤波器要保证阻带截止频率大于1000Hz，通带截止频率小于2500Hz； 而带通滤波器要保证阻带截止频率1大于500Hz，通带截止频率1小于1000Hz，通带截止频率2大于1000Hz，阻带截止频率2小于2500Hz。Bandpass模块的实际参数设置如图5.7.6所示。需要注意的是，Frequency units（频率单位）选择Hz，其值应与前文信号产生模块的抽样频率

图5.7.6Bandpass模块的参数设置

保持一致。设定好参数后，可以通过单击View Filter Response按钮观察滤波器模块的各种响应，如幅频响应，如图5.7.7所示。在此界面，可以将滤波器的分子、分母多项式系数保存到文本文档，以作他用。



图5.7.7Bandpass模块的幅频响应


仿真时间设置为1s，通过Scope模块观察输入、输出信号的波形，结果如图5.7.8和图5.7.9所示。


图5.7.8例5.7.3的输入信号





图5.7.9例5.7.3的3个滤波器的输出信号


5.8时域抽样和恢复
5.8.1信号的时域抽样

工程中，信号的时域抽样就是从连续时间信号x(t)中抽取出一系列离散样本值。抽样间隔一般是恒定的，记作Ts，那么抽样得到的离散样本值是


x(t)|t=nTs=x(nTs)Δx(n)
（5.28）


这里的x(n)即为离散时间信号。
例5.8.1设连续时间信号x(t)=Sa2(t)，抽样间隔Ts=0.2s，画出被抽信号和抽样信号的波形图。
解： MATLAB源代码如下

clc; clear all; close all;

t=-10:0.01:10; 

xt=sinc(t/pi).^2;%被抽信号

Ts = 0.2;%抽样间隔为0.2s

n = round(min(t)/Ts):round(max(t)/Ts); %离散时间信号的自变量

xn = sinc(n*Ts/pi).^2;%抽样信号

subplot(2,1,1); plot(t,xt); title('连续时间信号x(t)')%画图

subplot(2,1,2); stem(n*Ts,xn,'filled'); title('连续时间信号的抽样序列x(n)')

程序运行结果如图5.8.1所示。


图5.8.1连续时间信号的抽样


5.8.2抽样信号的频谱
为了便于理论分析，本书中一般将周期冲激串作为抽样脉冲，得到


xs(t)=x(t)·δT(t)=∑∞n=-∞x(nTs)δt-nTs
（5.29）


经过推导


Xsjω=1Ts∑∞
k=-∞X［j(ω-kωs)］
（5.30）


从而得出结论： 抽样信号的频谱是被抽信号频谱的周期延拓，延拓周期是抽样角频率ωs。若按照5.8.1节的方法进行抽样，离散时间序列x(n)的频谱是否仍旧是被抽信号频谱的周期延拓，延拓的周期是否仍旧是ωs呢？通过例5.8.2回答上述问题。
例5.8.2画出例5.8.1中被抽信号和抽样序列的频谱图。
解： MATLAB源代码如下

clc; clear all; close all;

%%%时域%%%%%%%%%

dt = 0.01;

t=-10:dt:10;%时域自变量

xt=sinc(t/pi).^2;%被抽样信号

Ts = 0.5;%抽样间隔

n = round(min(t)/Ts):round(max(t)/Ts); %离散时间信号的自变量

nTs = n*Ts;

xn = sinc(nTs/pi).^2;%抽样信号

%%%%频域%%%%%%%%

%被抽信号

w = -30:0.1:30; %频域自变量

［ww,tt］ = meshgrid(w,t);%变为矩阵形式

Xjw = dt*xt*exp(-j*tt.*ww);%矩阵-向量乘法计算被抽信号频谱

%抽样信号

［wws,tts］ = meshgrid(w,nTs);%变为矩阵形式

Xs_jw = Ts*xn*exp(-j*tts.*wws);

subplot(2,2,1); plot(t,xt);

xlabel('t'); title('连续时间信号x(t)')%画图

subplot(2,2,3); stem(nTs,xn,'.');

xlabel('t'); title('连续时间信号的抽样序列x(n)')

subplot(2,2,2); plot(w,Xjw);

xlabel('＼omega'); title('连续时间信号的频谱')

subplot(2,2,4); plot(w,Xs_jw);

xlabel('＼omega'); title('抽样序列的频谱')

xticks(［min(w) -2*pi/Ts 0 2*pi/Ts max(w)］)%在横坐标上显示ws

程序运行结果如图5.8.2所示。观察发现，抽样序列的频谱是被抽信号频谱的周期延拓，延拓的周期是抽样角频率ωs。只不过抽样序列频谱的主周期和被抽信号频谱相等，而不像式（5.30）那样有个系数1Ts。这是因为式（5.30）假设抽样脉冲是周期冲激串，而实际应用中并不是这样。


图5.8.2抽样前后信号的波形及频谱


5.8.3时域抽样定理
根据5.8.2节的介绍，抽样序列的频谱仍然是被抽信号频谱的周期延拓，且延拓的周期仍然是抽样角频率ωs，因此若要保证抽样序列保留了被抽信号的所有信息，仍然需要满足


ωs≥2ωm
（5.31）


其中,ωm是被抽信号的最高角频率。对应到抽样间隔Ts，需要满足


Ts≤πωm
（5.32）


例5.8.3利用例5.8.2，验证时域抽样定理。
解： 例5.8.2中，被抽信号x(t)=Sa2(t)，其傅里叶变换X(jω)=πΛ4(ω)，最高角频率ωm=2rad/s，故抽样间隔应小于π2s。为便于观察，保持被抽信号不变，不断增加抽样间隔，MATLAB源代码如下

clc; clear all; close all;

%%%被抽信号%%%%%%%%%

dt = 0.01;

t=-10:dt:10;%时域自变量

xt=sinc(t/pi).^2;%被抽样信号

w = -30:0.1:30; %频域自变量

［ww,tt］ = meshgrid(w,t);%变为矩阵形式

Xjw = dt*xt*exp(-j*tt.*ww);%矩阵-向量乘法计算被抽信号频谱

Ts = 0.3*pi;%抽样间隔，分别选取0.3π和0.6π

n = round(min(t)/Ts):round(max(t)/Ts); %离散时间信号的自变量

nTs = n*Ts;

xn = sinc(nTs/pi).^2;%抽样信号

［wws,tts］ = meshgrid(w,nTs);%变为矩阵形式

Xs_jw = Ts*xn*exp(-j*tts.*wws);

figure;

subplot(2,2,1); plot(t,xt);

xlabel('t'); title('连续时间信号x(t)')%画图

subplot(2,2,3); stem(nTs,xn,'filled');

xlabel('t'); title('抽样序列x(n)')

subplot(2,2,2); plot(w,Xjw);

xlabel('＼omega'); title('连续时间信号的频谱');xlim(［-20 20］)

subplot(2,2,4); plot(w,Xs_jw);

xlabel('＼omega');

hold on; plot(［-2 2］,［pi pi］,'r')

plot(［-2,-2］,［0,pi］,'r')

plot(［2,2］,［0,pi］,'r');

hold off

title(［'抽样序列的频谱,抽样间隔为',num2str(Ts/pi),'π'］);xlim(［-20 20］);

程序运行结果如图5.8.3和图5.8.4所示，动态结果请扫描二维码。为便于比较，图中用方框框出了［-2，2］的角频率范围。当抽样间隔不满足式（5.32）时，抽样序列的频谱将出现混叠现象。







图5.8.3抽样间隔为0.3π的结果





图5.8.4抽样间隔为0.6π的结果


5.8.4抽样信号的恢复
在满足时域抽样定理的前提下，理论教材的结论是可以将抽样信号经过模拟理想低通滤波器恢复出来。需要注意的是，工程中抽样序列是离散的，无法经过模拟滤波器，通常做法是将其通过数模转换器(DAC)。数模转换器主要采用零阶抽样保持或一阶抽样保持，示意图如图5.8.5和图5.8.6所示。


图5.8.5零阶抽样保持示意图




图5.8.6一阶抽样保持示意图


例5.8.4绘制抽样、恢复过程中信号的波形及频谱。
解： 对例5.8.2中的抽样信号分别进行零阶、一阶抽样保持。在MATLAB中，可以采用线性插值函数interp1完成一阶抽样保持，但没有直接实现零阶抽样保持的函数，需要自行编程实现。MATLAB源代码如下

clc; clear all; close all;

%%%时域%%%%%%%%%

dt = 0.01;

t=-10:dt:10;%被抽信号时域自变量

xt=sinc(t/pi).^2;%被抽样信号

Ts = 0.5;%抽样间隔

n = round(min(t)/Ts):round(max(t)/Ts); %离散时间信号的自变量

nTs = n*Ts;

xn = sinc(nTs/pi).^2;%抽样信号

nKeep = 50;%恢复时的保持点数

yt0 = repmat(xn,nKeep,1);%零阶抽样保持信号

yt0 = transpose(yt0(:));

nt = linspace(-10,10,length(yt0));%恢复信号的时域自变量

yt1 = interp1(nTs,xn,nt,'linear');%一阶抽样保持信号

%%%%频域%%%%%%%%

w = -30:0.1:30; %频域自变量

%被抽信号

［ww,tt］ = meshgrid(w,t);%变为矩阵形式

Xjw = dt*xt*exp(-j*tt.*ww);%矩阵-向量乘法计算被抽信号频谱

%抽样信号

［wws,tts］ = meshgrid(w,nTs);%变为矩阵形式

Xs_jw = Ts*xn*exp(-j*tts.*wws);%抽样序列频谱

%恢复信号

［wwr,ttr］ = meshgrid(w,nt);%变为矩阵形式

Yjw_0 = Ts/nKeep*yt0*exp(-j*ttr.*wwr);%零阶抽样保持信号的频谱

Yjw_1 = Ts/nKeep*yt1*exp(-j*ttr.*wwr);%一阶抽样保持信号的频谱

subplot(4,2,1); plot(t,xt);

xlabel('t'); title('被抽信号')%画图

subplot(4,2,3); stem(nTs,xn,'filled');

xlabel('t'); title('抽样序列')

subplot(4,2,5); plot(nt,yt0);

xlabel('t'); title('零阶抽样保持恢复信号')%画图

subplot(4,2,7); plot(nt,yt1);

xlabel('t'); title('一阶抽样保持恢复信号')

subplot(4,2,2); plot(w,Xjw);

xlabel('＼omega'); title('被抽信号的频谱')

subplot(4,2,4); plot(w,Xs_jw);

xlabel('＼omega'); title('抽样序列的频谱')

subplot(4,2,6); plot(w,Yjw_0);

xlabel('＼omega'); title('零阶抽样保持恢复信号的频谱')

subplot(4,2,8); plot(w,Yjw_1);

xlabel('＼omega'); title('一阶抽样保持恢复信号的频谱')

程序运行结果如图5.8.7所示。与被抽信号频谱相比，零阶抽样保持恢复信号的频谱在低频部分没有发生变化，而在高频部分多了一些成分，这是由于时域跳变引起的； 而一阶抽样保持恢复信号的频谱不仅低频部分保持较好，而且由于时域跳变较小，高频部分的变化也较小，相比而言恢复效果更好。


图5.8.7抽样、恢复信号