02 第2章 MATLAB在信号与系统中的应用.pptx

第1页：第2章MATLAB在信号与系统中的应用

在信号与线性系统分析中涉及的信号的傅里叶变换、拉普拉斯变换等内容是MATLAB在通信原理、数字信号处理、数字图像处理等课程中应用的基础。本章通过连续、离散信号的表示及可视化、信号的卷积和、信号的自相关函数、互相关函数、傅里叶变换、功率谱估计、拉普拉斯变换、系统的零极点分析等实例展示了在MATLAB中进行信号与系统知识的仿真的基本方法和基本原理。

第2页：本章知识架构

【本章知识架构】

第3页：典型连续信号的时域表示及可·实验目的

1．实验目的

(1)掌握信号的MATLAB表示方法。

(2)利用MATLAB实现信号的时域运算与时域变换。2．实验原理在MATLAB中通常用到两种方法来表示信号：一种是向量表示信号；另一种是用符号运算的方法来表示信号。用适当的MATLAB表示出信号后，就可以利用MATLAB的绘图命令绘出真正的信号波形。

第4页：例1·典型连续信号的时域表示及可

【例1】用MATLAB画出信号f(t)=ε(t+3)-2ε(t)的波形，如图2.3所示。 f=sym('heaviside(t+3)-2*Heaviside(t)') ezplot(f,[-5,4]),

hold on保留当前图形及坐标属性，以便下一条绘图语句继续在当前图形窗口上绘图

plot([0,0],[-1,1]), axis([-5,4,-1,1.1])

hold off恢复至默认模式，即plot语句在绘制下一个图形之前擦除原来的绘图及坐标属性

第5页：单位阶跃信号

在MATLAB中还提供一些图形函数，专门用于对plot命令所画出图形进行修饰。如坐标范围设定(axis命令)、加坐标轴名称(xlabel，ylabel命令)、给图形加标题(title命令)、绘网格(grid命令)、对图形进行文字标注(text命令)。图2.3单位阶跃信号

第6页：典型连续信号的时域表示及可视化（续）

在MATLAB中可以在同一个坐标中画出多幅图形。但plot命令在执行时首先将当前图像窗口清屏，用户看到的是最后一条命令绘制的图形。MATLAB提供了hold命令，将当前图形窗口的图形保留，利用多条plot命令绘制多幅图形。hold on命令：保留当前窗口的图形。hold off命令：解除hold on命令。

第7页：典型连续信号的时域表示及可·程序代码

注意

观察axis、title、text这3个语句的应用。

t=0:0.05:3; f=exp(-1.5*t); plot(t,f) axis([0,3,0,1.2]) title('单边指数信号') text(3.1,0.05,'t') (2)用符号函数sgn(t)来生成单位阶跃函数。ε(t)=1/2+(1/2)sgn(t)，而sgn(t)的表示可

调用MATLAB中的符号函数sign来实现：

t=-5:0.05:3; f=sign(t);

然后，用下面的命令就可以绘出符号函数的波形(执行结果如图2.5所示)：

plot(t,f) axis([-5,5,-1.1,1.1])

第8页：符号函数的波形

图2.5符号函数的波形

第9页：再用以下的命令表示单位阶跃信程序代码

再用以下的命令表示单位阶跃信号并绘出它的波形(执行结果如图2.6所示)。

ff=1/2+1/2*f; plot(t,ff),axis([-5,5,-1.1,1.1])图2.6用符号函数画出的单位阶跃信号

第10页：单位冲激信号

单位冲激信号

第11页：典型离散信号的时域表示及可·实验目的

1．实验目的

(1)掌握MATLAB语言的基本操作，学习基本的编程功能。

(2)熟悉连续信号经过采样前后的频谱变化，加深对采样定理的理解。2．实验原理一般来说，离散时间信号用f(k)表示，其中变量k为整数，代表离散的采样时间点。在MATLAB中，用一个向量f即可表示一个有限长度序列。但是这样的向量并没有包含其对应的时间序列信息。所以，要完整地表示离散信号需要用两个向量。

在MATLAB表示离散序列并将其可视化时，要注意以下几点：

• ①与连续时间信号不同，离散时间信号无法用符号运算来表示；
• ②由于在MATLAB中，矩阵的元素个数是有限的，因此MATLAB无法表示无限序列；
• ③在绘制离散信号波形时，要用专门的绘制离散数据的stem命令。

第12页：离散信号的图形显示使用程序代码

离散信号的图形显示使用stem指令。

编写MATLAB程序，产生下列典型脉冲序列。

(1)单位抽样序列。

(2)单位阶跃序列。

(3)正弦序列。

(4)实指数序列：

3x= (0.75)n。 (5)复指数序列：( 0.2 j0.7) 4x= e − +n。

4．程序代码

1)单位抽样序列

function dwxulie(k1,k2,k0) k=k1:k2; n=length(k); f=zeros(1,n); f(1,-k0-k1+1)=1 stem(k,f,'filled') axis([k1,k2,0,1.5]) title ('单位序列δ(k)')3．仿真思路

第13页：单位阶跃序列程序代码

2)单位阶跃序列

function dwjieyuexulie(k1,k2,k0) k=k1:-k0-1; kk=-k0:k2; n=length(k); nn=length(kk); u=zeros(1,n); uu=ones(1,nn); stem(kk,uu,'filled') hold on stem(k,u,'filled') hold off axis([k1,k2,0,1.5]) title ('单位阶跃序列')

第14页：正弦序列程序代码

3)正弦序列

function zhengxuanxulie() k=0:40; subplot(2,1,1) stem(k,cos(k*pi/8),'filled') title('cos(k*pi/8)') subplot(2,1,2) stem(k,cos(k*2),'filled') title('(cos(k*2))')

4)实指数、复指数序列

function zhishuxulie() n0=0;nf=10;ns=3; n3=n0:nf;x3=(0.75).^n3; %实指数序列 n4=n0:nf;x4=exp((-0.2+0.7j)*n4); %复指数冲序列 stem(n3,x3); %画出实指数序列的图 title('实指数序列') figure subplot(2,2,1),stem(n4,real(x4)); %复指数序列的实部 title('复指数序列实部') subplot(2,2,2),stem(n4,imag(x4)); %复指数序列的虚部 title('复指数序列虚部') subplot(2,2,3),stem(n4,abs(x4)); %复指数序列的幅度 title('复指数序列幅度') subplot(2,2,4),stem(n4,angle(x4)); %复指数序列的相位 title('复指数序列相位')

第15页：典型离散信号的时域表示及可·运行结果与分

5．运行结果与分析单位抽样序列如图2.9所示(函数调用方式为dwxulie(-5，5，3))。图2.9单位抽样序列

第16页：典型离散信号的时域表示及可视化（续）

单位阶跃序列如图2.10所示(函数调用方式为dwjieyuexulie(-5，5，3))。图2.10单位阶跃序列

第17页：正弦序列程序代码

正弦序列如图2.11所示(函数调用方式为zhengxuanxulie())。图2.11正弦序列

第18页：典型离散信号的时域表示及可视化（续）

实指数离散序列如图2.12所示(函数调用方式为zhishuxulie())。图2.12实指数离散函数

第19页：复指数序列实部

复指数序列实部、虚部、幅度、相位如图2.13所示。图2.13复指数序列实部、虚部、幅度、相位

第20页：典型离散信号的时域表示及可·思考题

6．思考题

(1)要得到ε(k−3)该调用哪个函数？如何调用？

(2)例中给出的两个正弦序列是周期序列吗？为什么？

第21页：信号的时域运算·实验目的

1．实验目的

• (1)了解信号常用的时域运算。
• (2)了解信号的时域运算在MATLAB中的实现方法。2．实验原理1)连续信号的时域运算与变换f1、f2是两个用符号表达式表示的连续信号，s、w、y均为运算或变换后得到的符号表达式。
• (1)相加

s=symadd(f1,f2)或s=f1+f2 ezplot(s)(2)相乘 w=symmul(f1,f2)或w=f1*f2 ezplot(w)

第22页：信号的时域运算·程序代码

(3)移位 y=subs(f,t,t-t0); ezplot(y)(4)反折 y=subs(f,t,-t); ezplot(y)(5)尺度变换 y=subs(f,t,a*t); ezplot(y)

对于以上命令，可在画图命令之后加入坐标轴的调整axis等命令，以使画出的图形更清晰、直观。2)离散序列的时域运算及变换对于离散序列来说，序列的相加、相乘是将两序列时间信号的值相加或相乘，平移等变换与连续信号的定义完全相同。但要注意的是，在MATLAB中，离散序列的相加、相乘需表示成两个序列的相加、相乘，因而参加运算的两序列必须有相同的维数。

第23页：信号的时域运算·仿真思路

3．仿真思路在序列求和的函数中，将要进行相加运算的两序列向量通过补零的方式成为同维数的两序列向量，因而在调用函数时，要进行相加运算的两序列向量维数可以不同。向量的反折可用MATLAB中的fliplr函数来实现。

第24页：信号的时域运算·程序代码

k=-3:3; f=2.^k; lsfz(f,k)

%注意比较两次出现波形的异同4)离散序列平移

function[f,k]=lsyw(ff,kk,k0) k=kk-k0;f=ff; subplot(2,1,1); stem(kk,ff,'filled') axis([min(min(kk),min(kk))-1, max(max(kk),max(k))+1,min(ff)-0.5,max(ff)+0.5]) ; title('原始序列'); subplot(2,1,2); stem(k,f,'filled') axis([min(min(kk),min(kk))-1, max(max(kk),max(k))+1,min(f)-0.5,max(f)+0.5]) title('平移后的序列');

第25页：信号的时域运算·运行结果与分析

5．运行结果与分析

第26页：信号的时域运算·运行结果图

图2.15序列相乘及其结果可视化图2.14序列相加及其结果可视化

第27页：序列反折及其可视化

序列反折及其可视化如图2.16所示。图2.16序列反折及其结果可视化

第28页：序列平移及其可视化

序列平移及其可视化如图2.17所示。为了更清楚地显示序列的变换，画图时对坐标轴的选取非常重要，这关系到结果可视化的直观性。注意分析在以上程序中，是如何选取坐标轴的范围的。图2.17序列平移及其结果可视化

第29页：

第30页：离散序列的卷积·实验目的

1．实验目的

(1)了解卷积和的概念。

(2)用conv函数计算卷积和。

(3)通过MATLAB仿真对卷积过程有感性认识。

第31页：

第32页：例5·离散序列的卷积

【例5】试用MATLAB计算如下算式序列的卷积和f(k)，绘出它们的时域波形。

f1(k)={1 2 1 0},f2(k)={1 1 1 1 1 0}

解：该问题可用dconv()函数来解决，实现这一过程的命令如下。

f1=[1 2 1 ]; k1=[-1 0 1]; f2=ones(1,5); k2=-2:2; [f,k]=dconv(f1,f2,k1,k2)4．运行结果与分析

离散序列的卷积和如图2.18所示。

第33页：离散序列的卷积（续）

由图可见，conv函数能正确计算两个离散序列的卷积和。本例的难点在于卷积和起止坐标的计算。只有正确计算坐标，画出来的图才是正确的。图2.18离散序列的卷积和

第34页：离散序列的卷积·思考题

5．思考题

(1)已知某LTI离散系统，其单位响应h(k)=ε (k)-ε (k−4)，求该系统在激励为f(k)=ε (k)− ε (k−3)时的零状态响应y(k)，并绘出其时域波形图。

(2)如果上题h(k)=ε (k)，重做上题。(提示：时间序列无限长时，必须将其进行截断，如只保留100个样值点)观察计算结果，所有计算样值均是真实的吗？

(3)编写连续系统卷积计算的子程序。

第35页：连续信号的卷积积分·实验目的

1．实验目的

(1)熟悉MATLAB软件平台。

(2)掌握MATLAB编程方法、常用语句和可视化绘图技术。

(3)编程实现常用信号及其运算MATLAB实现方法。2．实验原理连续信号是指自变量的取值范围是连续的，且对于一切自变量的取值，除了有若干个不连续点之外，信号都有确定的值与之对应。严格来说，MATLAB并不能处理连续信号，而是用等时间间隔点的样值来近似表示连续信号。当取样时间间隔足够小时，这些离散的样值就能较好地近似连续信号。由于MATLAB运算的特点，两个连续信号f1(t)、f2(t)的卷积f(t)=f1(t)*f2(t)，用MATLAB实现的过程如下。

(1)将连续信号f1(t)、f2(t)以时间间隔Δ进行取样，得离散序列f1(k)、f2(k)。

(2)构造f1(k)、f2(k)与相对应的时间向量k1和k2。

(3)调用conv()函数计算卷积积分f(t)的近似向量f(k)。

(4)构造f(k)对应的时间向量k。

第36页：

第37页：卷积子函数程序代码

2)卷积子函数sconv.m function [f,k]=sconv(f1,f2,k1,k2,p) %计算连续信号卷积积分f(t)=f1(t)*f2(t)

%f:卷积积分f(t)对应的非零值向量

%k:f(t)的对应时间向量

%f1,f2:f1(t),f2(t)的非零样值向量 %k1,k2:f1(t),f2(t)的对应时间向量

%p:取样时间间隔

f=conv(f1,f2); %计算序列f1,f2的卷积和f f=f*p; k0=k1(1)+k2(1); %计算序列f非零样值的起点位置 k3=length(f1)+length(f2)-2; %计算卷积和f的非零样值的宽度

k=k0:p:k3*p; %确定卷积和f非零样值的时间向量

subplot(2,2,1); plot(k1,f1); title('f1(t)'); xlabel('t'); ylabel('f1(t)'); subplot(2,2,2); plot(k2,f2); title('f2(t)'); xlabel('t'); ylabel('f2(t)'); subplot(2,2,3); plot(k,f); %画出卷积f(t)的波形 title('f(t)=f1(t)*f2(t))'); xlabel('t'); ylabel('f(t)');

第38页：连续信号的卷积积分·程序代码

h=get(gca,'position'); h(3)=2.5*h(3);

set(gca,'position',h) %将第3个子图的横坐标扩展为原来的2.5倍5．运行结果与分析信号及卷积结果图如图2.19所示。从上面结果能够看出f1(t)和f2(t)的波形以及它们的卷积波形，加深了理解卷积的数学运算过程。连续信号卷积与离散信号卷积过程从本质上来讲没有什么区别，都是调用conv函数来实现的。在连续信号中，要用plot函数作图，并且信号及卷积的结果的时间间隔都取得很小。

第39页：连续信号的卷积积分·运行结果图

图2.19信号及卷积结果图

第40页：

第41页：

第42页：循环卷积和线性卷积·程序代码

4．程序代码

1)卷积循环的主程序

n=[0:1:15]; m=[0:1:7]; N1=length(n); N2=length(m); xn=0.9.^n; %生成x(n) hn=ones(1,N2); %生成h(n) y1n=conv(xn,hn); %直接用函数conv计算线性卷积 y2n=circonv1(xn,hn,N1+N2-1); %用函数circonv1计算N1+N2-1点循环卷积 y3n=circonv1(xn,hn,N1); %用函数circonv1计算N1点循环卷积

第43页：循环卷积和线性卷积·程序代码

ny1=[0:1:length(y1n)-1]; ny2=[0:1:length(y3n)-1]; subplot(5,1,1); stem(n,xn);title('x(n)'); axis([0,25,0,1]); subplot(5,1,2); stem(m,hn);title('h(n)'); axis([0,15,0,1]); subplot(5,1,3); stem(ny1,y1n);title('x(n)与h(n)的线性卷积'); subplot(5,1,4); stem(ny1,y2n);title('x(n)与h(n)的N1+N2-1点循环卷积'); subplot(5,1,5); stem(ny2,y3n);title('x(n)与h(n)的N1点循环卷积'); axis([0,25,0,6]);

第44页：循环卷积circonv1程序程序代码

2)循环卷积circonv1程序 function y=circonv1(x1,x2,N) if length(x1)>N error('N must not be less than length of x1') end if length(x2)>N error('N must not be less than length of x2')

end %以上语句判断两序列的长度是否大于N

x1=[x1,zeros(1,N-length(x1))]; %填充x1(n),使其长度为N x2=[x2,zeros(1,N-length(x2))]; %填充x2(n),使其长度为N n=[0:1:N-1]; x2=x2(mod(-n,N)+1); %生成序列x2((-n))N H=zeros(N,N); for n=1:1:N H(n,:)=cirshiftd(x2,n-1,N); %该矩阵的n行x2((k-1-n))N end y=x1*H'; %计算卷积

第45页：循环移位cirshiftd程程序代码

3)循环移位cirshiftd程序 function y=cirshiftd(x,m,N) if length(x)>N error('length of x must be less than N'); end x=[x,zeros(1,N-length(x))]; %填充序列x(n)到指定长度 n=[0:1:N-1];

y=x(mod(n-m,N)+1); %利用取模实现循环移位5．运行结果与分析循环卷积和线性卷积如图2.20所示。第一幅和第二幅图是原始信号，第三幅图是直接用函数conv计算得到的线性卷积，是正确的结果；第四幅图是使用函数circonv1计算的(N1+N2−1)点循环卷积，从图中，可以看出它与第三幅图一致，这说明没有发生混叠；第五幅图是用函数circonv1计算的N1点的循环卷积，明显看出，它的结果的前面(N2−1)点与第三幅图不同，发生了混叠。由于采用快速傅里叶变换计算循环卷积要比直接计算线性卷积快很多，在实际中，计算两个序列的线性卷积都是通过把它化为循环卷积而实现的。用这个方法计算线性卷积，要通过合适的零值填充，使得结果不会发生混叠，从而得到正确的结果。

第46页：循环卷积和线性卷积·思考题

6．思考题

• (1)在什么情况下线性卷积与循环卷积的结果一样？为什么？
• (2)使用循环卷积的意义何在？图2.20循环卷积与线性卷积

第47页：自相关函数·实验目的

1．实验目的

(1)进一步了解自相关函数、互相关函数的基本理论。

(2)了解自相关函数的应用。

(3)学会如何对信号中加有白噪声的信号的自相关函数进行分析。2．实验原理

合理利用相关函数有利于检测和识别淹没在随机噪声中的周期信号。自相关定义为：

第48页：

第49页：自相关函数·仿真思路

3．仿真思路通过MATLAB仿真，查看一个加有白噪声的正弦信号的自相关函数以及白噪声的自相关函数，并因此了解数字通信系统的性能。MATLAB信号处理工具箱提供了计算随机信号相关函数的函数xcorr和计算协方差的函数xcov。函数xcorr用于计算随机序列的互相关函数和自相关函数。本例中主要使用xcorr函数，其调用格式为：

c=xcorr(x,y) c=xcorr(x,y,'option') c=xcorr(x,y,maxlags,'option') [c,lags]=xcorr(x,y,maxlags,'option')

第50页：自相关函数·程序代码

4．程序代码

N=1000; n=[0:N-1]; fs=500; t=n/fs; lag=100; x=sin(2*pi*8*t)+0.85*randn(1,length(t)); [c,lags]=xcorr(x,lag,'unbiased'); subplot(2,2,1); plot(t,x); xlabel('t'); ylabel('x(t)'); title('原始信号x'); grid; subplot(2,2,2); plot(lags/fs,c); xlabel('t'); ylabel('Rx(t)'); title('x的自相关函数'); grid; d=randn(1,length(x));

第51页：自相关函数·程序代码

ylabel('Rx(t)'); title('x的自相关函数'); grid; d=randn(1,length(x)); [c,lags]=xcorr(d,lag,'unbiased'); subplot(2,2,3); plot(t,d); xlabel('t'); ylabel('d(t)'); title('噪声信号d'); grid; subplot(2,2,4); plot(lags/fs,c); xlabel('t'); ylabel('Rd(t)'); title('d的自相关函数'); grid;

第52页：自相关函数·运行结果与分析

5．运行结果与分析不同信号的自相关函数如图2.21所示。由图2.21可见，含有周期成分和干扰噪声信号的自相关函数在t=0时具有最大值，且在t较大时仍具有明显的周期性，其频率和周期成分的周期相同；而不含周期成分的纯噪声信号在t=0时也具有最大值，但在稍大时明显至零。此性质可以用来识别随机信号中是否含有周期信号成分和它们的频率。从上述图形能简单地看出，自相关函数具有偶函数的性质。且白噪声信号在时间间隔大时，自相关函数趋向于0，即白噪声的均值，说明白噪声在时间间隔大时是无关的，这正表征了白噪声的随机性质。而周期性正弦信号在时间间隔大时其自相关并不趋向于0，是有关的。

第53页：自相关函数·思考题

6．思考题

• (1)自相关函数有什么用途？有什么特点？试根据图说明自相关函数的特点。
• (2) xcorr函数是用来计算互相关函数的，怎么样用它来计算自相关？
• (3)式(2-12)是如何得出的？图2.21不同信号的自相关函数

第54页：互相关函数·实验目的

1．实验目的

(1)进一步了解互相关函数的基本理论。

(2)了解互相关函数的应用。

(3)加强对MATLAB的进一步掌握。2．实验原理在信号处理过程中经常要研究两个信号的相似性，或者一个信号经过一段延迟后自身的相似性，以实现信号的检测、识别与提取等。相关函数是描述随机信号的重要统计量，已广泛应用于雷达、声呐、数字通信、地质学以及其他学科的工程领域中。

第55页：互相关函数·程序代码

4．程序代码

n=[1:50]; x= sin(pi/8*n+pi/4)+2*cos(pi/7*n); w=randn(1,length(n)); y=x+w; [c1,rxx]=xcorr(x); [c2,rxy]=xcorr(x,y); [c3,ryy]=xcorr(y); figure(1) subplot(221); plot(x); title('信号x波形图'); grid; subplot(222); plot(rxx,c1);

第56页：互相关函数·程序代码

title('信号x的自相关函数波形图');

grid; subplot(223); plot(y); title('信号y波形图'); grid; subplot(224); plot(ryy,c3);

title('信号y的自相关函数波形图');

grid; figure(2) subplot(111); plot(rxy,c2);

title('信号x和信号y的互相关函数波形图');

grid;5．运行结果与分析

信号的自相关函数如图2.22所示。

第57页：互相关函数·运行结果图

图2.22信号的自相关函数图2.23不同信号的互相关函数

第58页：程序的运行结果

程序的运行结果如图2.22、图2.23所示。显然，仅从信号y(n)的波形很难分辨出是否含有正余弦信号，而从它和x(n)的互相关函数以及自相关函数的波形中可以判别出信号y(n)含有正余弦分量。故相关函数可以进行噪声检测、信号中隐含的周期性检测等。6．思考题

(1)信号y的自相关函数有什么特点？信号x和y的互相关函数有什么特点？

(2)互相关函数经常用于什么地方？

第59页：周期信号的傅里叶级数·实验目的

1．实验目的

(1)学习周期信号傅里叶级数分解与综合方法。

(2)对周期信号进行频谱分析。2．实验原理周期信号是定义在(−∞，+ ∞)区间，按一定时间间隔(周期T)不断重复的信号。它可表示为：f(t) =f(t+mT) (2-19)式中：m为任意整数，T为周期，周期的倒数称为该信号的频率。

第60页：MATLAB实现上述分析过程的流程如下

MATLAB实现上述分析过程的流程如下。

• (1)编写子函数x=time_fun_x(t)，用符号表达式表示出周期信号在第一个周期内的符号表达式，并赋值给符号变量x。
• (2)编写子函数y=time_fun_e(t)，求出该周期信号在绘图区间内的信号样值，并赋值给变量y。
• (3)编写求解傅里叶系统及绘制合成波形图的通用函数CTFShchsym.m，这是一个计算连续时间周期信号的三角级数前7次展开系数，再用这7次系数合成原连续时间周期信号的程序。该函数的流程如下。

①调用函数time_fun_x(t)，获取周期信号的符号表达式。

②求出信号的傅里叶系数。

③求出各次谐波。

④绘制各次谐波叠加波形图。

⑤调用函数time_fun_e(t)，绘制原波形图。3．仿真思路

第61页：周期信号的傅里叶级数·程序代码

4．程序代码

1) CTFShchsym函数 function [A_sym,B_sym]=CTFShchsym

%采用符号计算求一个周期内连续时间函数f的三角级数展开系数，再用这些

%展开系数合成连续时间函数f傅里叶级数

%函数的输入输出都是数量值

% Nf=6谐波的阶数

% Nn输出数据的准确位数

% A_sym第一元素是直流项，其后元素依次是1,2,3…次谐波cos项展开系数

% B_sym第2,3,4…元素依次是1,2,3…次谐波sin项展开系数 % tao=1 tao/T=0.2 syms t n k x T=5;tao=0.2*T;a=0.5; if nargin<4;Nf=6;end

第62页：周期信号的傅里叶级数·程序代码

if nargin<5;Nn=32;end x=time_fun_x(t); AO=2*int(x,t,-a,T-a)/T; %求出三角函数展开系数AO As=int(2*x*cos(2*pi*n*t/T)/T,t,-a,T-a); %求出三角函数展开系数As Bs=int(2*x*sin(2*pi*n*t/T)/T,t,-a,T-a); %求出三角函数展开系数Bs A_sym(1)=double(vpa(AO,Nn)); %获取串数组AO所对应的ASC2码数值数组 for k=1:Nf A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); %获取串数组A所对应的ASC2码数值

数组

B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn)); %获取串数组B所对应的ASC2码数值数

组é

end if nargin==0

c=A_sym;disp(c) %输出为c三角级数展开系数：第一元素是直流项，其后元素依次

是1,2,3…次谐波cos项展开系数

d=A_sym;disp(d) %输出为d三角级数展开系数：第一元素是直流项,其后元素依次是

1,2,3…次谐波sin项展开系数

t=-8*a:0.01:T-a; f1=0.4/2+0.3742.*cos(2*pi*1*t/5)+0.*sin(2*pi*1*t/5); %基波

第63页：周期信号的傅里叶级数·程序代码

f2=0.3027.*cos(2*pi*2*t/5)+0.*sin(2*pi*2*t/5); %2次谐波 f3=0.2018.*cos(2*pi*3*t/5)+0.*sin(2*pi*3*t/5); %3次谐波 f4=0.0935.*cos(2*pi*4*t/5)+0.*sin(2*pi*4*t/5); %4次谐波 f5=-0.0624.*cos(2*pi*6*t/5)+0.*sin(2*pi*6*t/5); %6次谐波 f6=f1+f2; %基波+2次谐波 f7=f6+f3; %基波+2次谐波+3次谐波 f8=f7+f4+f5; %基波+2次谐波+3次谐波+4次谐波+6次谐波 subplot(2,2,1) plot(t,f1),hold on y=time_fun_e(t); %调用连续时间函数-周期矩形脉冲 length(y) length(t) plot(t,y,'r:');

title('周期矩形波的形成-基波')

axis([-4,4.5,-0.5,1.3]); subplot(2,2,2); plot(t,f6),hold on; y=time_fun_e(t); plot(t,y,'r:');

title('周期矩形波的形成-基波+2次谐波');

第64页：周期信号的傅里叶级数·程序代码

axis([-4,4.5,-0.5,1.3]); subplot(2,2,3); plot(t,f7),hold on y=time_fun_e(t); plot(t,y,'r:') title('基波+2次谐波+3次谐波') axis([-4,4.5,-0.5,1.3]) subplot(2,2,4) plot(t,f8),hold on y=time_fun_e(t); %调用连续时间函数-周期矩形脉冲 plot(t,y,'r:');

title('基波+2次谐波+3次谐波+4次谐波+6次谐波')

axis([-4,4.5,-0.5,1.3]) end

第65页：周期信号的傅里叶级数·程序代码

2) time_func_x函数 function x=time_fun_x(t)

%该函数是CTFShchsym.m的子函数,它由符号变量和表达式写成.

h=1; x1=sym('Heaviside(t-0.5)')*h; x=x1-sym('Heaviside(t-0.5)')*h;3) time_func_e函数 function y=time_fun_e(t)

%该函数是CTFShchsym.m的子函数，它由符号变量和表达式写成

a=0.5;T=5;h=1;tao=0.2*T; e1=1/2+1/2.*sign(t+tao/2); e2=1/2+1/2.*sign(t-tao/2); y=h.*(e1-e2); %连续时间函数-周期矩形脉冲

第66页：周期信号的傅里叶级数·运行结果与分析

5．运行结果与分析由图2.24可见，基波是构成矩形脉冲的基础部分；2次谐波与基波同相相加，使矩形脉冲的主体更突出，而拖尾部分得到衰减，更逼近矩形波；3次、4次、6次谐波更进一步修饰矩形脉冲的细节，使矩形脉冲拖尾更小，峰值更突出。由此可见，低频成分构成信号的主体，而高频成分修饰信号的细节。该性质对于绝大多数信号而言是成立的，该性质是语音信号及图像信号压缩理论的理论基础之一。

第67页：周期信号的傅里叶级数·思考题

6．思考题

• (1)上例的图2.24可以验证什么结论？对上例增加几次谐波又如何?
• (2)使用MATLAB绘出周期矩形脉冲信号的振幅频谱。图2.24周期矩形波的形成

第68页：

第69页：

第70页：

第71页：

第72页：信号的傅里叶变换·运行结果图

图2.25程序3的信号的幅度谱图2.26信号f(t)及其频谱

第73页：

第74页：信号的功率谱估计·实验目的

1．实验目的

(1)进一步了解功率谱估计原理。

(2)掌握经典法与现代法功率谱估计的各种计算方法。2．实验原理在通信系统中，往往需要研究具有某种统计特性的随机信号。由于随机信号是一类持续时间无限长、具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此，不能用确定信号的频谱计算方法去分析随机信号的频谱。然而，虽然随机信号的频谱不存在，但其相关函数却是确定的。如果随机信号是平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数，简称功率谱。功率谱反映了单位频带内随机信号的功率大小，它是频率的函数，其物理单位是W/Hz。随机信号功率谱估计就是根据随机信号的一个样本信号来对该随机信号的功率谱密度函数作出估计。对功率谱估计常用的方法有周期图法、AR模型法、MUSIC算法等。

第75页：信号的功率谱估计·仿真思路

3．仿真思路对高斯噪声中的信号x(t)=cos(2π×60×t)+cos(2π×63×t)进行功率谱估计，采样频率为fs=1000Hz，分别用周期图法、Yule-Walker法、Burg法、改进的协方差法与Music法估计序列的功率谱，且AR模型的阶数取30与50两种情况讨论。4．程序代码

fs=1000;%采样频率 t=0:1/fs:1;%含有噪声的序列 xn=cos(2*pi*60*t)+cos(2*pi*63*t)+0.1*randn(size(t));%计算序列DFT nfft=1024; window=boxcar(length(xn)); [Pxx,f]=periodogram(xn,window,nfft,fs); figure(1) plot(f,10*log10(Pxx)) grid on xlabel('Frequence/Hz'); ylabel('PSD/dB') title('Periodgram PDS Estimate')

第76页：信号的功率谱估计·程序代码

order1=30; order2=50; figure(2) pyulear(xn,order1,nfft,fs) figure(3) pyulear(xn,order2,nfft,fs) figure(4) pburg(xn,order1,nfft,fs) figure(5) pburg(xn,order2,nfft,fs) figure(6) pmcov(xn,order1,nfft,fs) figure(7) pmcov(xn,order2,nfft,fs)%Music法 p=50; figure(8) pmusic(xn,p,nfft,fs)

第77页：信号的功率谱估计·运行结果与分析

5．运行结果与分析图2.28周期图法功率谱估计图2.29 Yule-Walker法功率谱估计(阶数30)

第78页：信号的功率谱估计·运行结果图

图2.31 Burg法功率谱估计(阶数30)图2.33修正协方差法功率谱估计(阶数30)

第79页：Music法功率谱估计

图2.35 Music法功率谱估计

第80页：信号的功率谱估计（续）

由以上对比可知，周期图法的频率分辨率最低，而YuleWalker法比Burg法的分辨率要低，采用修正协方差法可以进一步提高分辨率，而且增加AR模型的阶数也可以提高分辨率。在几种方法中，Music法分辨率最高。6．思考题

(1)几种功率谱估计方法中，哪种分辨率最好，为什么？

(2)如何对周期图法进行改进？

第81页：不同加窗方式对功率谱估计·实验目的

1．实验目的

(1)研究不同的加窗方式对信号频谱估计的影响。

(2)加强对MATLAB的进一步掌握。2．实验原理连续非周期信号频谱的数值计算必须首先对信号进行时域采样，得到按时间离散化的信号，时域采样必须满足或近似满足采样定理。根据时域频域的对应关系，时域采样将导致采样信号频谱周期化。然而，为了使周期化后采样信号频谱便于计算机处理，还必须再将其进行频域离散化，方法是对该频谱进行频域采样。根据时域频域的对应关系，频域离散化将对应于时域信号的周期化。因此，对连续非周期信号频谱进行数值计算时，要确定如何截取信号的时间段、如何选择时域采样率，以及在时间段上对信号进行截取的方式。截取信号的时间段长度决定了时域周期化的周期，对应于频域采样的频率间隔即频率分辨率；时域采样率决定了频域周期化的周期即频谱数值计算的范围；而在某时间段上对信号进行截取的方式，即不同窗函数的应用，决定了信号频谱估计的精度和有效范围。

第82页：不同加窗方式对功率谱估计·仿真思路

3．仿真思路假定要分析连续时间非周期信号f(t)在频率范围[0，fm]内的频谱，且要求分析的频率分辨率为FHz，则首先应根据信号频率范围确定采样率，再根据所要求的频率分辨率确定截取时间长度，从而计算出所需计算FFT的序列长度即点数，最后，根据信号时域波形特征选择使用不同的窗函数。

第83页：不同加窗方式对功率谱估计·程序代码

4．程序代码

clear all; fs=200; F=1; T=1/fs; L=1/F; N=floor(fs/F)+1; t=0:T:L; freq=0:F:fs; ft=(1.2*sin(2*pi*60*t)+0.8*sin(2*pi*80*t))'; ft_rectwin=rectwin(N).*ft./sqrt(sum(abs(rectwin(N).^2))./N); ft_hamming=hamming(N).*ft./sqrt(sum(abs(rectwin(N).^2))./N); ft_hann=hann(N).*ft./sqrt(sum(abs(rectwin(N).^2))./N); Fw_rectwin=T.*fft(ft_rectwin,N); Fw_hamming=T.*fft(ft_hamming,N);

第84页：不同加窗方式对功率谱估计·程序代码

Fw_hann=T.*fft(ft_hann,N); figure(1) subplot(2,2,1);plot(t,ft); title('原信号'); subplot(2,2,2);plot(t,ft_rectwin); title('矩形窗'); subplot(2,2,3);plot(t,ft_hamming); title('汉明窗'); subplot(2,2,4);plot(t,ft_hann); title('汉宁窗'); figure(2) subplot(3,1,1);plot(freq,10*log10(abs(Fw_rectwin))); title('矩形窗频谱');ylabel('幅度dB'); axis([0 200-50 0]);grid on; subplot(3,1,2);plot(freq,10*log10(abs(Fw_hamming))); title('汉明窗频谱');ylabel('幅度dB'); axis([0 200-50 0]);grid on; subplot(3,1,3);plot(freq,10*log10(abs(Fw_hann))); title('汉宁窗频谱');ylabel('幅度dB'); axis([0 200-50 0]);grid on;

第85页：不同加窗方式对功率谱估计·运行结果图

图2.36加窗后的时域波形图2.37加窗后的信号幅度谱

第86页：不同加窗方式对功率谱估计·运行结果与分析

5．运行结果与分析由仿真结果可以看出，矩形窗的优点是频谱的主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄露；汉明窗与汉宁窗都是余弦窗，只是加权系数不同，汉明窗加权的系数能使旁瓣达到更小，但旁瓣衰减速度比汉宁窗慢；汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄露观点出发，汉宁窗优于矩形窗和汉明窗，但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨率下降。

第87页：不同加窗方式对功率谱估计·思考题

6．思考题

(1)信号功率谱估计时为什么要加窗？不同的加窗方式对信号功率谱的估计有何影响？

(2)频域周期化的周期即频谱数值计算的范围是由什么决定的？信号功率谱估计的精度和范围是由什么决定的？

第88页：基于Welch法的功率谱估·实验目的

1．实验目的进一步了解功率谱估计方式。

2．实验原理Welch法是一种使用最为广泛的经典功率谱估计方法，该方法在周期图法的基础上主要进行了两处改进：其一是将采样数据进行分段，并允许各段数据有交叠的部分；其二是每段数据加窗处理时，窗函数不局限于矩形窗，也可以用汉宁窗或海明窗等，从而可以改善谱失真现象。

第89页：基于Welch法的功率谱估·仿真思路

3．仿真思路

(1)先将N个数据分为L段，每段M个数据。

(2)选择适当的窗函数w(n)，对每段数据依次加权，然后确定每段的周期图。

(3)对分段周期图进行平均，得到功率谱估计。4．程序代码

MATLAB代码示例如下。

clear; Fs=1000; n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); nfft=1024; window=boxcar(100); %矩形窗 window1=hamming(100); %汉明窗 window2=blackman(100); %blackman窗 noverlap=20; range='half'; %频率间隔为[0 Fs/2]，只计算一半的频率

第90页：基于Welch法的功率谱估·程序代码

[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range); [Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range); [Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range); plot_Pxx=10*log10(Pxx); plot_Pxx1=10*log10(Pxx1); plot_Pxx2=10*log10(Pxx2); figure(1) plot(f,plot_Pxx); pause; figure(2) plot(f,plot_Pxx1); pause; figure(3) plot(f,plot_Pxx2);

第91页：基于Welch法的功率谱估·运行结果与分

5．运行结果与分析由运行结果可见，Welch法更适合于用FFT进行计算。主要提出两方面的修正：其一是选择适当的窗函数ω(n)，并在周期图计算前直接加进去，这样加窗函数的优点是无论什么样的窗函数均可使谱估计非负。其二是在分段时，可使各段之间有重叠，这样将会使方差减小(当N与M一定时)。利用MATLAB实现Welch算法，从而完成功率谱估计，不仅实现方便，形象直观，还有利于分析窗函数对功率谱估计的影响。总之，要针对不同的信号和不同的处理目的来选择合适的窗函数，这样才能得到良好的效果。

第92页：基于Welch法的功率谱估·运行结果图

图2.38加矩形窗的功率谱估计图2.39加汉明窗的功率谱估计

第93页：基于Welch法的功率谱估·思考题

6．思考题

• (1)经典谱估计算法有什么缺点？与经典谱估计方法相比，Welch算法做了什么改进？
• (2) Welch算法中对窗函数的选择有什么要求？图2.40加Blackman窗的功率谱估计

第94页：

第95页：拉普拉斯变换及其曲面图·程序代码

4．程序代码

1)单位阶跃信号的拉普拉斯变换

x1=-0.2:0.03:0.2; y1=-0.2:0.03:0.2; [x,y]=meshgrid(x1,y1); s=x+i*y;

fs=abs(1./s); %计算拉普拉斯变换在复平面上的取样点

mesh(x,y,fs); %绘制的拉普拉斯变换曲面图 surf(x,y,fs);

title('单位阶跃信号拉普拉斯变换曲面图');

colormap(hsv); axis([-0.2,0.2,-0.2,0.2,0.60]); rotate3d;

第96页：单边正弦信号程序代码

2)单边正弦信号f(t)=sin(t)ε (t)的拉氏变换曲面图

clf; a=-0.5:0.08:0.5; b=-1.99:0.08:1.99; [a,b]=meshgrid(a,b); d=ones(size(a)); c=a+j*b; c=c.*c; c=c+d; c=1./c; c=abs(c); mesh(a,b,c); surf(a,b,c); axis([-0.5,0.5,-2,2,0,15]);

title('单边正弦信号拉普拉斯变换曲面图')；

colormap(hsv);

第97页：拉普拉斯变换及其曲面图·运行结果与分析

5．运行结果与分析图2.41单位阶跃信号拉普拉斯变换曲面图图2.42单边正弦信号拉普拉斯变换曲面图

第98页：拉普拉斯变换及其曲面图·思考题

6．思考题求下列信号的拉普拉斯变换，并绘制拉普拉斯变换在S平面的三维曲面图。

(1)f(t)=cos(2t)ε(t) (2)f(t)=sin(πt)[ε(t)−ε(t−1)]

第99页：离散系统的零极点分析·实验目的

1．实验目的

(1)通过MATLAB仿真简单的离散时间系统，研究其时域、频域特性。

(2)加深对离散系统的冲激响应、频率响应分析和零极点分布概念的理解。

第100页：离散系统的零极点分析·仿真思路

3．仿真思路

(1)计算离散时间系统的频率响应函数：[h,f] = freqz(b,a,n,fs)。

(2)系统的有理分式形式转化为零极点增益形式的函数：[z,p,k] = tf2zp(b,a)。

(3)绘制离散系统零极点图函数：zplane(b,a)、zplane(z,p,k)。

• (4)例子，设离散时间系统对应的输入输出差分方程为：y(n) − 0.4y(n−1) − 0.5y(n− 2) = 0.2x(n) + 0.1x(n−1)
• ①求系统的频率响应。
• ②画出系统的零极点图

第101页：离散系统的零极点分析·程序代码

4．程序代码

1)求系统的频率响应的程序

fs=1000; %取样频率 b=[0.2 0.1]; %系统函数分子的系数 a=[1-0.4-0.5]; %系统函数分母的系数 [h,f]=freqz(b,a,256,fs); %计算离散时间系统的频率响应 mag=abs(h); %取出幅度 ph=angle(h); %取出相位 ph=ph*180/pi; subplot(2,1,1),plot(f,mag); %画图 title('幅频特性'); xlabel('频率(Hz)'); ylabel('幅度'); subplot(2,1,2),plot(f,ph); title('相频特性'); xlabel('频率(Hz)'); ylabel('相位');

第102页：画出零极点的程序程序代码

2)画出零极点的程序

b=[0.2 0.1]; a=[1-0.4-0.5]; [z,p,k]=tf2zp(b,a);

zplane(z,p); %画出以零点矢量z和极点矢量p以及增益k描述的离散时间的零极点图

title('零极点分布图')

或者：

b=[0.2 0.1]; a=[1-0.4-0.5];

zplane(b,a); %画出以矢量b和a描述的离散时间系统的零极点图

title('零极点分布图')

第103页：离散系统的零极点分析·运行结果图

图2.43离散时间系统的幅频特性与相频特性图2.44离散时间系统的零极点分布图

第104页：线性系统的根轨迹分析·实验目的

1．实验目的分析系统根轨迹与系统响应之间的关系。2．实验原理闭环控制系统的稳定性和性能指标由闭环系统的极点在复平面的位置决定，因此分析和设计系统时确定出闭环极点的位置非常有意义。要求出闭环的极点，就要求解高阶代数方程，三阶以上的代数方程求解是非常困难的，而且，每当参数变化时，还要对代数方程反复求解，因此为了解决此问题，人们一直在寻求不用解方程就能确定某参数变化时极点的位置的方法。于是提出了根轨迹法，此方法可以由开环传递函数的零极点求出闭环的极点。所谓根轨迹是指系统的开环的传递函数中的某个参数(如开环增益K)从零变到无穷时，闭环的特征根在S平面上移动的轨迹；当变化的参数是开环增益时，所对应的根轨迹为常规根轨迹；当变化的参数为开环中其他的参数时对应的轨迹为广义的根轨迹；一般而言，当变化系统为正反馈时，对应的根轨迹为零度根轨迹，而负反馈系统的根轨迹为180º根轨迹。

第105页：

第106页：阶跃响应曲线为了进一步分析系程序代码

2)阶跃响应曲线为了进一步分析系统的稳定性，分别绘制K=2.5、K=6、K=7时，系统闭环单位阶跃响应曲线。

num=1,den=[1 3 2 0]; %传递函数的分子系数向量和分母系数向量

num1=2.5*num; %当K=2.5时，传递函数分子系数向量 sys1=tf(num1,den); %建立传递函数模型 sys11=feedback(sys1,1); %求闭环传递函数 figure(1); %设置第一个窗形窗口 step(sys11); %求系统sys11的单位阶跃响应曲线 title('k=2.5单位阶跃响应') %添加图名称 num2=6*num; %当K=6时,传递函数分子系数向量 sys2=tf(num2,den); %建立传递函数模型 sys22=feedback(sys2,1); %求闭环传递函数 figure(2); %设置第二个窗形窗口 step(sys22); %求系统sys22的单位阶跃响应曲线 title('k=6单位阶跃响应') %添加图名称 num3=7*num; %当K=7时,传递函数分子系数向量 sys3=tf(num3,den); %建立传递函数模型 sys33=feedback(sys3,1); %求闭环传递函数 figure(3); %设置第三个窗形窗口 step(sys33); %求系统sys33的单位阶跃响应曲线

title('k=7单位阶跃响应¦') %添加图名称

第107页：线性系统的根轨迹分析·运行结果与分析

5．运行结果与分析分析图2.46，判断系统稳定性。程序执行后，可在图形窗口根轨迹中显示当选择根轨迹上某点时，其相应增益由变量K记录，相应极点由poles记录，可见当参数K(0→6)变动时，根轨迹均在S平面纵坐标左侧，对应系统是稳定的，一旦根轨迹穿越纵坐标到达右侧，对应K> 6，系统就不稳定，如图2.47、图2.48所示。

第108页：线性系统的根轨迹分析·运行结果图

图2.46系统的根轨迹图图2.47K=2.5时的单位阶跃响应

第109页：线性系统的根轨迹分析·运行结果图

图2.48K=6时的单位阶跃响应图2.49K=7时的单位阶跃响应

第110页：线性系统的根轨迹分析（续）

由上面的3幅图可见，当K=2.5时，系统单位阶跃响应曲线为衰减振荡；当K=6时，系统单位阶跃响应曲线为等幅振荡，对应的系统为临界稳定；当K=7时系统单位阶跃响应曲线为发散振荡，对应的系统不稳定。这完全验证了以上的运算结论。6．思考题系统是否稳定与系统的根轨迹有什么关系？

第111页：本章小结

本章小结本章通过16个实例展示了信号与线性系统分析课程中核心知识点用MATLAB仿真实现的基本方法及基本知识。本章既根据信号与线性系统分析课程的知识框架给出了对应的仿真实例，也是后面MATLAB在通信原理、数字信号处理、数字图像处理等课程应用的基础。

第112页：结束页

Thanks书山有路勤为径，学海无涯苦作舟。好好学习，天天向上