信号与系统 第二版 重大 例题程序 第8章

例8.1用Matlab命令绘制连续时间信号。Matlab程序如下：t=0:0.01:10;%设定时间向量f=5*exp(-0.5*t).*cos(2*t);%产生信号f(t)plot(t,f)%绘制f(t)xlabel('t'),ylabel('f(t)')%设定坐标轴名程序运行结果如图：图8.2.1例8.1信号波形例8.2用Matlab命令绘制冲击函数信号。Matlab程序如下：t=-5:0.01:5;%设定时间向量y=0*(t>=-5&t<0)+1*(t==1)+0*(t>0&t<=5);%根据冲击函数性质给出冲击函数表达式plot(t,y)%绘制冲击函数axis([-5,5,-1,1.5]);%设置坐标范围程序运行结果如图：图8.2.2例8.2信号波形例8.3用Matlab命令绘制阶跃函数信号。在Matlab中，产生阶跃函数有多种方法：用符号函数实现；用stepfun()实现；用比较表达式实现。（1）用符号函数实现的Matlab程序如下：t=-200:200;%设定时间向量f1=sign(t);%产生符号函数f2=1/2+f1/2;%产生单位阶跃信号subplot(1,2,1),plot(t,f1)%绘制符号函数axis([-200200-1.11.1])%设置坐标范围xlabel('t'),ylabel('sgn(t)')%设置坐标轴名title('符号函数')%设置图名subplot(1,2,2),plot(t,f2)%绘制单位阶跃信号axis([-200200-0.11.1])%设置坐标范围title('单位阶跃信号')%设置图名xlabel('t'),ylabel('u(t)')%设置坐标轴名程序运行结果如图：图8.2.3例8.3用符号函数实现的信号波形（2）stepfun（T,T0）：T为时间向量，T0为阶跃时刻，当T<T0时返回0，当T>T0时返回1，且返回的向量与T具有相同的长度。用stepfun()实现的Matlab程序如下：t=linspace(-200,200,40000);%设定时间向量y=stepfun(t,1);%设定1为阶跃时刻plot(t,y);%绘制阶跃函数xlabel('t'),ylabel('x'),%设置坐标轴名axis([-200200-0.11.1])%设置坐标范围程序运行结果如图：图8.2.4例8.3用stepfun()实现实现的信号波形（3）用比较表达式实现的Matlab程序如下：t=linspace(-200,200,40000);%设定时间向量x=t>0;%设置比较时刻plot(t,x)%绘制阶跃函数xlabel('t'),ylabel('x'),%设置坐标轴名axis([-200200-0.11.1])%设置坐标范围程序运行结果如图：图8.2.5例8.3用比较表达式实现实现的信号波形例8.4用Matlab求解例2.6中的LTI系统的全响应。解：先求零输入响应，再求零状态响应，两部分加起来就是系统的全响应。1）Matlab程序如下：eq1='D2y+5*Dy+6*y=0';%设定零输入条件下的微分方程ic1='y(0)=2,Dy(0)=-1';%输入初始状态yzi=dsolve(eq1,ic1);%求解微分方程，得到零输入响应yzi=simplify(yzi)%化简零输入响应eq2='D2y+5*Dy+6*y=2*exp(-t)';%设定给定输入条件下的微分方程ic2='y(0)=0,Dy(0)=0';%设定初始状态为0yzs=dsolve(eq2,ic2);%求解微分方程，得到零状态响应yzs=simplify(yzs)%化简零状态响应y=yzi+yzs%全响应运行程序，得到：yzi=exp(-3*t)*(5*exp(t)-3)yzs=exp(-3*t)*(exp(t)-1)^2y=exp(-3*t)*(5*exp(t)-3)+exp(-3*t)*(exp(t)-1)^2全响应结果经化解整理后与例2.6第1）求得的全解一致。2）Matlab程序如下：eq1='D2y+5*Dy+6*y=0';%设定零输入条件下的微分方程ic1='y(0)=1,Dy(0)=0';%输入初始状态yzi=dsolve(eq1,ic1);%求解微分方程，得到零输入响应yzi=simplify(yzi)%化简零输入响应eq2='D2y+5*Dy+6*y=exp(-2*t)';%设定给定输入条件下的微分方程ic2='y(0)=0,Dy(0)=0';%设定初始状态为0yzs=dsolve(eq2,ic2);%求解微分方程，得到零状态响应yzs=simplify(yzs)%化简零状态响应y=yzi+yzs%全响应运行程序，得到：yzi=exp(-3*t)*(3*exp(t)-2)yzs=exp(-3*t)*(t*exp(t)-exp(t)+1)y=exp(-3*t)*(3*exp(t)-2)+exp(-3*t)*(t*exp(t)-exp(t)+1)全响应结果经化解整理后与例2.6第2）求得的全解一致。例8.5用Matlab求解第2章例2.9中的LTI系统的零输入响应、零状态响应和全响应。解：此处采用求解零状态响应函数initial()、系统的时域响应函数lsim()进行编程作图求解。Matlab程序如下：a=[132];%设定系数向量ab=[26];%设定系数向量bsys=tf(b,a);%得到系统的传递函数模型[A,B,C,D]=tf2ss(b,a);%由传递函数模型求得状态空间模型sys1=ss(A,B,C,D);%状态空间模型y0=[20];%输入系统的初始条件t=0:0.01:5;f=heaviside(t);%系统输入y1=initial(sys1,y0,t);%零输入响应subplot(2,2,1);plot(t,y1);%绘制多子图xlabel('t');ylabel('y1')title('零输入响应')y2=lsim(sys1,f,t);%零状态响应subplot(2,2,2);plot(t,y2);xlabel('t');ylabel('y2')title('零状态响应')y3=y1+y2;%全响应subplot(2,1,2);plot(t,y3);xlabel('t');ylabel('y3')title('全响应')程序运行结果如下图：图8.2.6例8.5各种响应曲线例8.6用Matlab求解第3章例3.1中的周期为T，脉冲幅度为1的矩形波信号的傅里叶级数展开。解：对矩形波信号进行傅里叶级数展开最关键的是求得傅里叶级数的系数、，Matlab编程需要制定周期为固定值，故这里假设周期为4，脉冲宽度为2。程序如下：T=4;%设定矩形脉冲信号的周期w=2*pi/T;%基频F=@(t)abs(t)<=0;a0=quad(F,-2,2)/T;%求a0N=50;%设定最大谐波数an=zeros(1,N);bn=zeros(1,N);%初始化和fork=1:Nan(k)=quad(@rectcos,-2,2,[],[],k,w)*2/T;%计算anbn(k)=quad(@rectsin,-2,2,[],[],k,w)*2/T;%计算bnend%绘制周期矩形脉冲信号t=-4:0.01:4;%设定时间向量x=pulstran(t,-6:4:6,'rectpuls',2);%生成周期矩形脉冲信号subplot(3,2,1);plot(t,x);%绘制周期矩形脉冲信号xlabel('t');ylabel('Amplitude');axis([-4,4,-0.2,1.2]);gridon;title('Rectangularwave');%有限项级数逼近k_max=[110304050];%最大谐波数向量form=1:length(k_max)f=a0;fork=1:k_max(m)f=f+an(k)*cos(k*w*t)+bn(k)*sin(k*w*t);endsubplot(3,2,m+1);plot(t,f);gridon;xlabel('t');ylabel('PartialSum');axis([-4,4,-1,1.2]);title(['MaxHarmonics=',num2str(k_max(m))]);end仿真结果如下图：图8.2.7例8.6有限项傅里叶级数和例8.7用Matlab求解第3章例3.17中均匀冲击序列的傅里叶变换。解：从图3.6.1中可得均匀冲激序列的波形图和频谱图。根据欧拉公式，可将展开成三角函数形式的傅里叶级数。级数中的每一项都是连续函数，这里将用傅里叶级数的有限项来近似单位冲激序列，观察有限项级数信号叠加后逐渐演变成离散冲激的过程。Matlab实现代码如下：T=1;%设周期为1w1=2*pi/T;t=-1.2*T:0.01:1.2*T;%设定时间起始值的间隔N=5;%取傅里叶级数的前5项近似单位冲激序列xN=ones(1,length(t))/T;%傅里叶级数首项表达式plot(t,xN,'m');%作首项图xlabel('t');%标识X轴ylabel('x_N(t)');%标识Y轴holdon%叠绘color=['m','r','g','b'];%设定曲线颜色集forn=1:N-1xn=cos(n*w1*t)/T*2;xN=xn+xN;%求傅里叶级数的后续项和plot(t,xN,color(mod(n,length(color))+1));%作前n项和图endh=plot(t,xN,'k');%作前5项和图set(h,'linewidth',3);%设定前5项和的曲线线宽为3axis([-1.2*T,1.2*T,min(xN)*1.1,max(xN)*1.1]);%设定坐标轴范围(a)N=5(b)N=20图8.2.8例8.7有限项傅里叶级数合成单位冲激序列例8.8求的傅里叶变化，并绘制幅度谱。解：傅里叶变换指令较简单，值得注意的一点的需要先定义符号变量，阶跃函数在Matlab中的表示为Heaviside。Matlab实现傅里叶变换代码如下：symst;%定义符号变量tf=sym('heaviside(t+1)')-sym('heaviside(t-1)');%输入符号函数fF=fourier(f);%求f的傅立叶变换Fezplot(abs(F));%绘制幅度谱图8.2.9例8.8幅度谱例8.9用Matlab求解第3章例3.21中LTI系统的频率响应。解：本题编辑程序的重点在于准确给出LTI系统的常微分方程的系数，然后使用freqs()指令即可通过图形展示该系统的幅频特性和相频特性。程序如下：b=[143];a=[12];%设定分子、分母多项式系数向量w=-20:0.05:20;%设定频率向量H=freqs(b,a,w);%求频率响应mag=abs(H);%取频率响应的模phase=angle(H)*180/pi;%取频率响应的相位，并转换成角度subplot(2,1,1),plot(w,mag)%以线性坐标绘制幅频特性曲线xlabel('Frequency(rad/s)'),ylabel('Magnitude'),gridsubplot(2,1,2),plot(w,phase)%以线性坐标绘制相频特性曲线xlabel('Frequency(rad/s)'),ylabel('Phase(degrees)'),grid图8.2.10例8.9频率响应曲线例8.10用Matlab求例4.4、例4.5信号的单边拉普拉斯变换。解：例4.4Matlab程序如下：symsta;%定义符号变量t，af=exp(-a*t).*heaviside(t);%输入符号函数F=laplace(f);%求f的拉普拉斯变换FF运行结果为:F=1/(a+s)例4.5Matlab程序如下：symsta;%定义符号变量t，af=1+exp(-a*t).*heaviside(t);%输入符号函数fF=laplace(f);%求f的拉普拉斯变换FF运行结果为：F=1/(a+s)+1/s注：对比Matlab运行结果与例4.4的计算结果，Matlab没有对收敛域进行区分，只是给出了积分收敛时的结果。例4.5的计算结果与Matlab运行结果一致。另外，这里没有给出a的具体值，所以不能使用ezplot()绘图。例8.11用Matlab求例4.24象函数的原函数。解：求原函数即求拉普拉斯逆变换，采用ilaplace()函数实现。Matlab程序如下：symss;%定义符号变量sL=(s+6)/(s*(s+1)*(s+2));%输入符号函数Lf=ilaplace(L);%求L的拉普拉斯逆变换ffezplot(f)%原函数作图运行结果为：f=2*exp(-2*t)-5*exp(-t)+3图8.2.11例8.11原函数曲线注：对比例4.24的计算结果与Matlab的仿真结果，Matlab仿真结果中没有阶跃函数，因为laplace()和ilaplace()函数均有单边变换，默认取值大于0。另外，由于该题原函数中各参数确定，可以用ezplot()函数作图直观表示原函数。例8.12用Matlab求例4.29LTI系统的零输入响应、零状态响应和全响应。解：由例4.29可知，对系统的微分方程两边取拉普拉斯变换可得：整理得并代入初始状态得：则可以用Matlab求其全响应、零输入响应、零状态响应。程序如下：symstsYzs=(s+4)/(s^2+5*s+6);yzs=ilaplace(Yzs)ft=heaviside(t);F=laplace(ft);Yzf=F*2*(s+3)/(s^2+5*s+6);yzf=ilaplace(Yzf)yt=simplify(yzs+yzf)运行结果为：yzs=2*exp(-2*t)-exp(-3*t)yzf=1-exp(-2*t)yt=exp(-2*t)-exp(-3*t)+1例8.13用Matlab求例5.5离散系统差分方程，，，激励为，求系统的完全响应。程序如下：a=[6,-5,1];b=1;n=0:20;x=10*cos(0.5*pi*n);%生成激励信号x（n）wi=[-10-a(2)*(1/2)-a(1)*(-5/4),0];%定义初始状态，wi=f(2)-a(2)*y(1)-a(1)*y(2),从题目中给出%的初始条件可以得知f(2)=-10，y(2)=-5/4[y,wf]=filter(b,a,x,wi);%滤波得到完全响应plot(n,y);xlabel('n');ylabel('y')运行结果如下：图8.3.5例8.13系统完全响应图例8.14用Matlab求例5.8离散系统差分方程的单位序列响应，并与例5.8理论值做对比。程序如下：

k=0:10;a=[1,-3,3,-1];b=[1];h=impz(b,a,k);subplot(2,1,1);stem(k,h);title('单位冲击响应近似值');gridon;hk=1+3/2.*k+1/2.*k.^2;%在时为1subplot(2,1,2);stem(k,hk);title('单位冲击响应理论值');gridon;运行结果如下：图8.3.6例8.14离散时间系统单位序列响应对比两个图可知，用求得的单位序列响应近似值与理论值基本相符。例8.15用Matlab求例5.10离散系统差分方程的单位序列响应和单位阶跃响应，并与理论值、做对比。程序如下：k=0:10;a=[1,-5/6,1/6];b=[1];h=impz(b,a,k);g=stepz(b,a,k);subplot(2,2,1)stem(k,h);title('单位序列响应近似值');gridon;subplot(2,2,2)stem(k,g);title('单位阶跃响应近似值');gridon;hk=3*(1/2).^k-2*(1/3).^k;subplot(2,2,3)stem(k,hk);title('单位序列响应理论值')gk=3-3*(1/2).^k+(1/3).^k;subplot(2,2,4)stem(k,gk);title('单位阶跃响应理论值')图8.3.7例8.15离散时间系统单位序列响应及单位阶跃响应对比上图可知，用求得的单位序列响应近似值与理论值基本相符，用求得的单位阶跃响应近似值与理论值基本相符。例8.16用Matlab求例5.11中，的卷积和。程序如下：n=[0:20];h1=exp(-n);h2=(n>=0);y=conv(single(h1),single(h2));%conv指令要求两个向量是一种精度，故需要进行转换；卷积和%的长度length(y)=length(h1)+length(h2)-1。n=0:2*length(n)-2;%stem绘图要求两个参数长度一致，故进行n的转换stem(n,y);xlabel('n');ylabel('y');图8.3.8例8.16卷积和图形例8.17用Matlab求例6.1中单位序列、例6.3中因果序列的z变换。（1）在Matlab中，kroneckerDelta(n,0)代表单位序列函数，参考程序如下：f=sym('kroneckerDelta(n,0)');F=ztrans(f)运行结果为：F=1程序运行结果与例6.1理论计算结果一致。（2）参考程序如下：f=sym('a^n');F=ztrans(f)运行结果为：F=-z/(a-z)程序运行结果与例6.3理论计算结果一致。例8.18用Matlab求例6.20中象函数的原序列。程序如下：F=sym('(z^2+z)/(z-1)^2');f=iztrans(F)运行结果如下：f=2*n+1注：利用iztrans指令求得的是收敛域时的原序列，需要在程序运行结果后加单位阶跃函数，结果为f=（2*n+1），与理论结果一致。例8.19用Matlab求例6.26中二阶离散系统差分方程，，的零输入响应、零状态响应和全响应。解：首先对差分方程左右两边取单边z变换，得整理得：接下来可以利用Matlab求其完全响应、零状态响应和零输入响应，程序如下：symsnzf(z)=2.^n;F(z)=ztrans(f(z))Ys=((5-6*z^-1)-6)/(1-5*z^-1+6*z^-2);ys=iztrans(Ys)Yf=z.^-1*F(z)/(1-5*z.^-1+6*z.^-2);yf=iztrans(Yf)yn=simplify(ys+yf)运行结果为：F(z)=z/(z-2)ys=8*2^n-9*3^nyf=3*3^n-4*2^n-2^n*(n-1)yn=5*2^n-2^n*n-6*3^n即得系统的零输入响应为：零输入响应为：系统的全响应为：例8.20用Matlab判断例6.28中离散系统的稳定性，系统的差分方程为程序如下：a=[1];b=[1/3,1/3,1/3];zplane(b,a);运行结果如下：图8.3.9例8.19离散系统的零极点图从图中可以看出，该系统有两个零点，两个极点，极点为0，在单位圆内，故系统稳定，与理论计算结果一致。例8.21用Matlab求例6.28中离散系统的系统函数，并判断稳定性，系统的差分方程为程序如下：a=[1,0.1,-0.2];b=[1,1];[r,p,k]=tf2zp(b,a)zplane(b,a);运行结果为：r=-1p=-0.50000.4000k=1图8.3.10例8.20离散系统的零极点图由运行结果知：该系统的零点为-1和0，极点为-0.5和0.4，系数为1，故得该系统的系统函数为：，由于极点都在单位圆内，故系统稳定，与例6.29理论计算结果一致。例8.22分别设计20阶低通模拟chebyshevI型滤波器，通带内的最大衰减为0.3dB；20阶低通模拟chebyshevII型滤波器，阻带内的最小衰减为45dB，并给出其频率特性图。程序如下：[z1,p1,k1]=cheb1ap(20,0.3);[num1,den1]=zp2tf(z1,p1,k1);%将零极点形式转换为系统函数形式[z2,p2,k2]=cheb2ap(20,45);[num2,den2]=zp2tf(z2,p2,k2);figure(1)freqs(num1,den1)figure(2)freqs(num2,den2)运行结果如下图：图8.4.3低通模拟chebyshevI型滤波器图8.4.4低通模拟chebyshevII型滤波器例8.23设计数字信号处理系统，采样频率，希望设计成butterworth型高通数字滤波器，通带内允许的最大衰减为0.4dB，阻带内的最小衰减为40dB，通带上线临界频率为400Hz，阻带下限频率为300Hz。解：设计数字高通滤波器的流程如下图8.4.5：图8.4.5数字高通滤波器设计流程Matlab程序如下：wp=400*2*pi;ws=300*2*pi;rp=0.4;rs=40;fs=2000;%转换成模拟滤波器[N,Wc]=buttord(wp,ws,rp,rs,'s');%计算阶数和截止频率[Z,P,K]=buttap(N);%产生模拟低通滤波器[A,B,C,D]=zp2ss(Z,P,K);%零极点形式转换为状态方程形式[AT,BT,CT,DT]=lp2hp(A,B,C,D,Wc);%低通向高通的转换[num1,den1]=ss2tf(AT,BT,CT,DT);%状态方程形式向传递函数形式转换[num2,den2]=bilinear(num1,den1,fs);%双线性法得到的滤波器[H,W]=freqz(num2,den2);%离散系统频率响应F=W*fs/2/pi;semilogy(F,abs(H));grid%表示y坐标轴是对数坐标系xlabel('频率（Hz）');ylabel('幅值');运行结果如下：图8.4.6例8.23数字高通滤波器仿真结果例8.24以中国移动的上行频率为例（下行频段设计类似），设计数字带通滤波器，设抽样频率为，使用椭圆型滤波器。要求：①通带范围为885-909Hz，在通带边缘频率出的衰减不大于0.5dB；②在780Hz以下和1000Hz以上衰减不小于15dB。Matlab程序如下：fs=2500;ws=[885,909]/(fs/2);wp=[780,1000]/(fs/2);rp=0.5;rs=15;[N,wc]=ellipord(wp,ws,rp,rs,'s');%计算阶数和截止频率[num,den]=ellip(N,rp,rs,wc);%设计椭圆滤波器[H,W]=freqz(num,den);F=W*fs/2/pi;plot(F,20*log10(abs(H)));gridxlabel('频率（Hz）');ylabel('幅值');axis([02000-500])运行结果如下：图8.4.7例8.24带通滤波器仿真结果例8.25用矩形窗设计线性相位FIR低通滤波器，通带的截止频率为0.2，滤波器单位脉冲响应的长度为21。Matlab程序如下：M=21;wc=0.2*pi;r=(M-1)/2nr=-r:r;hdn=sin(wc*nr)/pi./nr;ifrem(M,2)~=0,hdn(r+1)=wc/pi;endw=boxcar(M);hn=hdn.*w';%相乘的时域冲激响应hw=fft(hn,512);w1=2*[0:511]/512;%频谱特性subplot(211),stem(1:M,hn,'o');xlabel('n'),ylabel('h(n)');subplot(212),plot(w1,20*log10(abs(hw)));xlabel('w/pi'),ylabel('幅值（dB）');axis([02-9010]);运行结果如下图：图8.4.8例8.25FIR低通滤波器仿真结果例8.26设计一个多通带滤波器，归一化的通带是：[0,0.2]、[0.4,0.6]、[0.8,1]。注意高频端为通带，故滤波器的阶数应该为偶数，取N=30。Matlab程序如下：wn=[0.2,0.4,0.6,0.8];N=30;b=fir1(N,wn,'dc-1');freqz(b)figure(2)stem(b,'.');xlabel('n'),ylabel('h(n)');运行结果如下：图8.4.9例8.26FIR多通带滤波器频谱特性图8.4.10例8.26单位冲激响应例8.27用fir2设计【例8.26】中的滤波器。Matlab程序如下：f=0:0.01:1m(1:21)=1;m(22:41)=0;m(42:61)=1;m(62:81)=0;m(82:101)=1;plot(f,m,'k');holdonN=30;b=fir2(N,f,m);[h,f1]=freqz(b);f1=f1./pi;plot(f1,abs(h));legend('理想滤波器','设计滤波器');xlabel('归一化频率'),ylabel('幅值');图8.4.11例8.27滤波器结果例8.28已知三角脉冲信号为，试用Matlab实现该信号经冲激脉冲抽样后得到的抽样信号及其频谱。解：设采用的抽样间隔为时，程序如下：Ts=1;dt=0.01;t1=-4:dt:4;xt=(4-abs(t1));subplot(221)plot(t1,xt),gridonaxis([-4.24.2-0.14.1]);xlabel('t'),title('三角脉冲信号')N=500;k=-N:N;W=pi*k/(N*dt);Xw=dt*xt*exp(-j*t1'*W);subplot(222)plot(W,abs(Xw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('三角脉冲信号频谱')t2=-4:Ts:4;xst=(4-abs(t2));subplot(223)plot(t1,xt),holdonstem(t2,xst),gridonaxis([-4.24.2-0.14.1]),xlabel('t'),title('经过冲激抽样后的信号'),holdoffXsw=Ts*xst*exp(-j*t2'*W);subplot(224)plot(W,abs(Xsw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('抽样后信号的频谱')运行结果如下：图8.4.13例8.28三角脉冲信号的冲激抽样例8.29试用三角脉冲信号来验证抽样定理。解：例8.27中三角脉冲信号的截止频率为，故奈奎斯特间隔，在例8.27中可以通过修改的值来得到不同的结果。取时和取时结果如下：Ts=1;dt=0.01;t1=-4:dt:4;xt=(4-abs(t1));subplot(421)plot(t1,xt),gridonaxis([-4.24.2-0.14.1]);xlabel('t'),title('三角脉冲信号')N=500;k=-N:N;W=pi*k/(N*dt);Xw=dt*xt*exp(-j*t1'*W);subplot(422)plot(W,abs(Xw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('三角脉冲信号频谱')t2=-4:Ts:4;xst=(4-abs(t2));subplot(423)plot(t1,xt),holdonstem(t2,xst),gridonaxis([-4.24.2-0.14.1]),xlabel('t'),title('经过冲激抽样后的信号'),holdoffXsw=Ts*(4-abs(t2))*exp(-j*t2'*W);subplot(424)plot(W,abs(Xsw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('抽样后信号的频谱')t3=-4:2:4;subplot(425)plot(t1,xt),holdonstem(t3,(4-abs(t3))),gridonaxis([-4.24.2-0.14.1]),xlabel('t'),title('t=2经过冲激抽样后的信号'),holdoffXsw=2*(4-abs(t3))*exp(-j*t3'*W);subplot(426)plot(W,abs(Xsw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('t=2抽样后信号的频谱')t4=-4:3:4;subplot(427)plot(t1,xt),holdonstem(t4,(4-abs(t4))),gridonaxis([-4.24.2-0.14.1]),xlabel('t'),title('t=3经过冲激抽样后的信号'),holdoffXsw=3*(4-abs(t4))*exp(-j*t4'*W);subplot(428)plot(W,abs(Xsw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('t=3抽样后信号的频谱')运行结果如下：图8.4.14例8.29三角脉冲信号的冲激抽样（Ts取不同的值）观察图8.4.14可以发现，抽样时间间隔大于奈奎斯特间隔时，信号的频谱会产生十分严重的混叠现象。例8.30对于三角脉冲信号，假设截止频率为，抽样间隔为，试用Matlab恢复抽样信号，并计算恢复后的信号与原信号的绝对误差。解：去低通滤波器的截止频率为，则Matlab程序如下：wm=pi/2;wc=1.1*wm;Ts=1;n=-50:50;nTs=n*Ts;xs=(4-abs(nTs));t=-4:0.01:4;xt=xs*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));t1=-4:0.01:4;x1=4-abs(t1);subplot(311)stem(nTs,xs),holdonplot(t1,x1),gridonaxis([-4.14.1-0.14.1])title('Ts=1时的抽样信号')holdoffsubplot(312)plot(t,xt),gridonaxis([-4.14.1-0.14.1])title('恢复后所得的三角脉冲信号')error=abs(xt-x1);subplot(313)plot(t,error),gridontitle('恢复信号与原信号的绝对误差')运行结果如下：图8.4.15例8.30抽样信号的恢复例8.31对信号、进行幅度调制。程序如下：fs=100;fc=15;t=0:0.01:2;x=sin([pi*t',2*pi*t']);y=ammod(x,fc,fs);z=amdemod(y,fc,fs);subplot(211)plot(t,x(:,1),'-',t,z(:,1),'o',t,y(:,1),'r')title('sin(pi*t)')subplot(212)plot(t,x(:,2),'--',t,z(:,2),'*',t,y(:,2),'r')title('sin(2*pi*t)')运行结果如下：图8.4.16例8.31幅度调制信号例8.32对信号、进行频率调制和解调。程序如下：fs=2000;fc=750;t=[0:fs]'/fs;s1=sin(2*pi*300*t)+2*sin(2*pi*600*t);s2=sin(2*pi*150*t)+2*sin(2*pi*900*t);x=[s1,s2];dev=50;y=fmmod(x,fc,fs,dev);z=fmdemod(y,fc,fs,dev);subplot(211)plot(t,x(:,1),'-',t,z(:,1),'o',t,y(:,1),'r')title('s1')subplot(212)plot(t,x(:,2),'--',t,z(:,2),'*',t,y(:,2),'r')title('s2')运行结果如下：图8.4.17例8.32频率调制信号例8.33对模拟信号进行相位调制，经信道叠加白噪声，解调并绘制原始信号与解调后信号波形。程序如下：fs=100;fc=10;t=[0:2*fs+1]'/fs;x=sin(2*pi*t)+sin(4*pi*t);phasedev=pi/2;y=pmmod(x,fc,fs,phasedev);y=awgn(y,10,'measured',103)%叠加噪声z=pmdemod(y,fc,fs,phasedev);figure;plot(t,x,'k-',t,z,'r-');legend('originalsignal','recoveredsignal')运行结果如下：图8.4.18例8.33相位调制解调信号例8.34给一锯齿波信号添加白噪声，并绘制出原信号和含噪声信号的波形。程序如下：t=0:0.1:10;x=sawtooth(t);y=awgn(x,10,'measured');plot(t,x,'*',t,y,'r')legend('originalsignail','signalwithawgn')运行结果如下：图8.4.18例8.34添加高斯白噪声的信号与原信号波形例8.35利用Simulink仿真实现对随机二进制信号（即随机出现的矩形脉冲函数）x(t)用cos200t进行调制，叠加随机噪声，再进行相干解调。观察各部分信号的情况及噪声对解调后的信号产生的影响。系统仿真模型图如下图8.4.19所示。图8.4.19例8.35随机矩形脉冲信号调制-解调Simulink仿真模型例8.36读取电脑中一段已经存在一段wav格式的语音信号语音一。程序如下：clc%清空当前工作区域[data,fs]=audioread('C:\Users\G\Documents\录音\ly1.wav');%读取数据left=data(:,1);right=data(:,2);%分离出左/右声道数据time=(1/fs)*length(left);%划分时间间隔t=linspace(0,time,length(left));%给定时间范围及步进plot(t,data)%做出频谱图像;soundsc(data,fs)%音频播放xlabel('time(sec)');ylabel('relativesignalstrength');%定义坐标轴表示运行结果可播放音频，频谱如下图8.4.22：图8.4.22音频信号频谱图例8.37对语音信号一和语音信号二进行预处理。程序如下：[y1,fs]=audioread('C:\Users\G\Documents\录音\ly1.wav');%读取语音一信号[y2,fs]=audioread('C:\Users\G\Documents\录音\ly2.wav');%读取语音二信号L1=length(y1);%测定语音一信号长度L2=length(y2);%测定语音二信号长度a1=y1(:,1).*hamming(L1);%加窗预处理a2=y2(:,1).*hamming(L2);%加窗预处理L1=length(a1);%测定语音一信号长度L2=length(a2);%测定语音二信号长度%采样信号的时域显示figure(1);subplot(211);plot(a1);title('语音一载波信号时域波形');subplot(212);plot(a2);title('语音二调幅信号时域波形');%傅里叶频谱绘制F1=fft(a1,L1);F2=fft(a2,L2);AF1=abs(F1);AF2=abs(F2);figure(2);subplot(211);plot(AF1);title('语音一载波信号幅频特性显示');subplot(212);plot(AF2);title('语音二调幅信号幅频特性显示');figure(3);freqz(F1);title('语音一载波信号FFT频谱显示');figure(4);freqz(F2);title('语音二载波信号FFT频谱显示');图8.4.23信号预处理之后时域波形图8.4.24信号预处理之后频域波形图8.4.25载波信号FFT频谱显示例8.38对例8.37中完成预处理的信号获取初始位置并进行信号合成。程序如下：%获取语音一信号的开始位置fori=1:L1-4g(i)=a1(i).*a1(i+1).*a1(i+2).*a1(i+3).*a1(i+4);%认为连续4个幅值不为0的信号即为开始ifg(i)~=0break;elsei=i+1;endendI=i;%获取语音二信号开始位置forj=1:L2-4m(j)=a2(j).*a2(j+1).*a2(j+2).*a2(j+3).*a2(j+4);ifm(j)~=0break;elsej=j+1;endendJ=j;%语音二信号hilbert变换H=hilbert(a2);figure(5);plot(abs(H));title('语音二信号包络显示');%信号对齐,语音二包络调制语音一振幅max1=max(I,J);fork=1:5*max1N(k)=a1(i).*H(j);i=i+1;j=j+1;end%N=N';N=N/(max(abs(N))*1.05);audiowrite('HC.wav',N,fs);%存储合成语音soundsc(imag(N),fs)%播放合成信号figure(6);plot(imag(N));title('合成信号时域显示');pause(1);FN=fft(N);figure(7);freqz(FN);title('合成声音信号FFT显示');figure(8);plot(abs(FN));title('合成声音信号的幅频特性');图8.4.26语音二信号包络图合成信号的时域显示结果、幅频特性、快速傅里叶变换结果分别如图8.4.27-8.4.29所示，该合成信号是以语音一信号的特性和语音二信号的幅度变化的，由其快速傅里叶变换的结果更证实了这一点，其幅频特性与语音一信号的幅频特性更接近。图8.4.27合成语音信号的时域波形图8.4.28合成语音信号的幅频特性图8.4.29合成语音信号快速傅里叶变换结果4、重建的趣味语音functionct2%定义常数FL=80;%帧长WL=240;%窗长P=10;%预测系数个数hw=hamming(WL);%汉明窗[s,Fs]=audioread('C:\Users\G\Documents\录音\ly1.wav');%载入语音sload('mtlb.mat');s=mtlb;s=buffer(s,FL);%分割信号x成不重叠的长度为FL的数据帧[~,FN]=size(s);%计算帧数%预测和重建滤波器exc=zeros(FL,FN);%激励信号（预测误差）s_rec=exc;%重建语音zi_pre=zeros(P,1);%预测滤波器的状态zi_rec=zeros(P,1);%重建滤波器的状态%合成滤波器exc_syn=exc;%合成的激励信号（脉冲串）s_syn=exc;%合成语音last_syn=0;%存储上一个帧的最后一个脉冲在下一帧结束的位置zi_syn=zeros(P,1);