通信原理MATLAB仿真教程 中国通信学会普通高等教育“十二五”规划教材立项项目 教学课件 ppt 作者 赵鸿图

信号与系统的MATLAB仿真*信号与系统的MATLAB仿真 4.1信号时域分析 4.2线性系统时域分析 4.3信号与系统的频域分析4.1信号时域分析 4.1.1常用信号的MATLAB表示 4.1.2信号运算 连续信号是指自变量的取值范围是连续的，且对于一切自变量的取值，除了有若干个不连续点之外，信号都有确定的值与之对应。严格来说，MATLAB并不能处理连续信号，而是用等时间间隔点的样值来近似表示连续信号。当取样时间间隔足够小时，这些离散的样值就能较好地近似连续信号。 因此，在MATLAB中，无论是连续时间信号还是离散时间信号，MATLAB都是用一个数字序列来表示信号，这个数字序列在MATLAB中叫做向量(vector)。通常的情况下，需要与时间变量相对应。 在MATLAB可视化绘图中，对于以t为自变量的连续信号，在绘图时统一用plot函数；而对n为自变量的离散序列，在绘图时统一用stem函数。一、几种常用信号的表示1．单位冲击函数、序列单位冲激函数、单位冲激序列 t=-5:0.01:5;y=(t==0);subplot(1,2,1);plot(t,y,'r');n=-5:5;x=(n==0);subplot(1,2,2);stem(n,x);2．单位阶跃函数、序列(一）几种常用信号的表示单位阶跃函数、单位阶跃序列 t=-5:0.01:5; y=(t>=0); subplot(1,2,1); plot(t,y,'r') n=-5:5;x=(n>=0);subplot(1,2,2);stem(n,x);3．实指数序列直接实现：n=[ns:nf];x=a.^n;4．复指数序列直接实现：n=[ns:nf];x=exp((sigema+jw)*n);5．正(余)弦序列直接实现：n=[ns:nf];x=cos(w*n+sita);生成上述三种信号 t=-5:0.01:5; subplot(2,2,1); a=2 y1=2.^t plot(t,y1,'r') subplot(2,2,2); a=2; theat=pi/3; y2=sin(2*pi*t+theat) plot(t,y2) subplot(2,1,2); w=4; y3=exp((a+j*w)*t); plot(t,y3,‘y')工具箱中的信号产生函数 函数名 功能 函数名 功能 sawtooth 产生锯齿波或三角波信号 pulstran 产生冲激串 square 产生方波信号 rectpule 产生非周期的方波信号 sinc 产生sinc函数波形 tripuls 产生非周期的三角波信号 chirp 产生调频余弦信号 diric 产生Dirichlet或周期sinc函数 gauspuls 产生高斯正弦脉冲信号 gmonopuls 产生高斯单脉冲信号 vco 电压控制振荡器1、sawtooth函数产生锯齿波或三角波例：产生f=50Hz的锯齿波、三角波Fs=10000;%采样频率t=0:1/Fs:0.1;%采样间隔1/Fsf=50;%50Hzx1=sawtooth(2*pi*50*t,0);x2=sawtooth(2*pi*50*t,1);x3=sawtooth(2*pi*50*t,0.5);subplot(311);plot(t,x1);subplot(312);plot(t,x2);subplot(313);plot(t,x3);2、square函数产生矩形波例：产生50Hz占空比分别为20％和50％的矩形波。Fs=10000;%采样频率t=0:1/Fs:0.1;%采样间隔1/Fsf=50;%50Hzx1=square(2*pi*50*t,20);x2=square(2*pi*50*t,50);subplot(211);plot(t,x1);subplot(212);plot(t,x2);3、sinc函数产生sinc波形或sin(πt)/(πt)波形例：t=linspace(-10,+10,200);x=sinc(t);plot(t,x);4.1.2信号运算1信号的相加与相乘y(n)=x1(n)+x2(n)MATLAB实现：y=x1+x2y(n)=x1(n)×x2(n)MATLAB实现：y=x1+x2;y=x1.*x2 function[f,k]=sigadd(f1,k1,f2,k2) %实现序列f1,f2的相加,相减,相乘,可据实际需要作选择 %f1,k1;f2,k2是参加运算的序列向量及其时间向量 %f,k作为返回的和(差,积)序列及其时间向量 %将f1,f2转换成等长序列s1,s2 k=min(min(k1),min(k2)):max(max(k1),max(k2)); s1=zeros(1,length(k));s2=s1;%初始化序列 s1(find((k>=min(k1))&(k<=max(k1))==1))=f1; s2(find((k>=min(k2))&(k<=max(k2))==1))=f2; f=s1+s2;%序列相加 %f=s1-s2;%序列相减 %f=s1.*s2;%序列相乘 stem(k,f,'fill')； axis([(min(min(k1),min(k2))-1),(max(max(k1),max(k2))+1),(min(f)-0.5),(max(f)+0.5)])2序列移位与周期延拓运算序列移位：y(n)=x(n-m)。MATLAB实现：y=x;ny=nx+m序列圆周移位：y(n)=x((n))M，MATLAB实现：ny=nxs:nxf；y=x(mod(ny,M)+1)3序列翻褶序列翻褶：y(n)=x(-n)。MATLAB可实现:y=fliplr(x)f(t+2)，f(-t)、f（（n-2））8f1=f;t1=t-2;subplot(2,2,2);plot(t1,f1)axis([-5502]);f2=fliplr(f);subplot(2,2,3);plot(t,f2).4两序列的卷积运算两序列卷积运算：MATLAB实现：y=conv(x1,x2)。序列x1(n)和x2(n)必须长度有限。连续卷积和离散卷积的关系： 所以，可以用离散卷积和CONV（）求连续卷积，只需足够小以及在卷积和的基础上乘以两序列的相关运算两序列相关运算：MATLAB实现：y=xcorr(x1,x2)。序列累加的数学描述为：MATLAB实现：y=cumsum(x)序列累加运算信号的能量和功率1.信号能量数字定义：MATLAB实现:E=sum(x.*conj(x));或E=sum(abs(x).^2);数字定义：2.信号功率MATLAB实现:P=sum(x.*conj(x))/N;或E=sum(abs(x).^2)/N;4.2线性系统时域分析 4.2.1线性时(移)不变系统特性 4.2.2线性时(移)不变系统表示 4.2.3线性时(移)不变系统的时域响应4.2.1线性时(移)不变系统特性 假设系统在输入信号x1(t)作用时的响应信号为y1(t)，在输入信号x2(t)作用时的响应信号为y2(t)，给定两个常数a和b，如果当输入信号为x(t)时系统的响应信号为y(t)，且满足 x(t)=x1(t)+x2(t)（4-15） y(t)=y1(t)+y2(t)（4-16） 则该系统具有叠加性（Additivity）。如果满足 x(t)=ax1(t)（4-17） y(t)=ay1(t)（4-18） 则该系统具有齐次性（Homogeneity）。一个系统如果是线性系统的话，那么这个系统必须同时具有叠加性和齐次性。4.2.2线性时(移)不变系统表示方法1．常系数线性微分/差分方程2．系统传递函数：tf3．零－极点增益模型zpk连续系统：连续系统：离散系统：离散系统：4．极点留数模型离散系统：连续系统：5．二次分式模型sos连续系统：离散系统：6．状态空间模型ss连续系统：离散系统：系统模型的转换函数在MATLAB中，用sos、ss、tf、zp分别表示二次分式模型、状态空间模型、传递函数模型和零－极点增益模型。其中sos表示二次分式，g为比例系数，sos为L×6的矩阵，即（4－15）1．ss2tf函数 格式：[num,den]=ss2tf(A,B,C,D,iu)功能：将指定输入量iu的线性系统(A,B,C,D)转换为传递函数模型[num,den]。2．zp2tf函数格式：[num,den]=zp2tf(z,p,k)功能：将给定系统的零－极点增益模型转换为传递函数模型，z、p、k分别为零点列向量、极点列向量和增益系数。线性系统模型的变换函数 函数名 功能说明 函数名 功能说明 ss2tf 状态空间模型转换为传递函数模型 zp2tf 零－极点增益模型转换为传递函数模型 ss2zp 状态空间模型转换为零－极点增益模型 zp2ss 零－极点增益模型转换为状态空间模型 ss2sos 状态空间模型转换为二次分式模型 zp2sos 零－极点增益模型转换为二次分式模型 tf2ss 传递函数模型转换为状态空间模型 sos2tf 二次分式模型转换为传递函数模型 tf2zp 传递函数模型转换为零－极点增益模型 sos2zp 二次分式模型转换为零－极点增益模型 tf2sos 传递函数模型转换为二次分式模型 sos2ss 二次分式模型转换为状态空间模型4.2.3线性时(移)不变系统的时域响应 一、连续系统的时域响应及MATLAB实现 二、离散时间系统的时域响应及MATLAB实现连续系统的时域响应及MATLAB实现 功能：求连续系统的单位冲激响应。 impulse（）函数有如下几种调用格式： impulse(b,a)：该调用格式以默认方式由向量a和b表示的连续系统在指定时间范围内的冲激响应h(t)的时域波形图，并能求出指定时间范围内冲激响应的数值解。 impulse(b,a,t)：绘出系统在0～t时间范围内冲激响应的时域波形。 impulse(b,a,t1:p:t2)：绘出在t1~t2时间范围内，且以时间间隔p均匀取样的冲激响应波形。 y=impulse(b,a,t1:p:t2)：不绘出波形，而是求出系统冲激响应的数值解：1、impulse函数连续系统的时域响应2、step函数功能：求连续系统的单位阶跃响应，与impulse（）函数一样，也有四种调用格式。3、initial函数功能：求连续系统的零输入响应。4、lsim函数功能：求任意输入的连续系统进行仿真。【例4-19】已知二阶系统y″(t)+y(t)=f(t);y′(0+)=-1，y(0+)=0，f(t)=cos2π,求系统的单位阶跃响应、零输入响应及全响应。 clear; b=[1];a=[101];%系统的传递函数参数 [ABCD]=tf2ss(b,a); sys=ss(A,B,C,D);%系统的状态空间 t=0:0.1:30;%响应的时间区间 f=cos(t);%系统的输入 zi=[-1;0];%系统的初始状态 y1=step(sys,t);%系统的单位阶跃响应 y2=initial(sys,zi,t);%系统的零输入响应 y3=lsim(sys,f,t);%系统的零状态响应 y4=lsim(sys,f,t,zi);%系统的全响应 subplot(141);plot(t,y1); xlabel('时间t');title('系统的单位阶跃响应'); line([0,30],[0,0]); subplot(142);plot(t,y2); xlabel('时间t');title('系统的零输入响应'); line([0,30],[0,0]); subplot(143);plot(t,y3); xlabel('时间t');title('系统的零状态响应'); line([0,30],[0,0]); subplot(144);plot(t,y4); xlabel('时间t');title('系统的全响应'); line([0,30],[0,0]);*1、impz函数求离散系统（数字滤波器）的单位冲击响应。离散时间系统响应·impz(b,a)以默认方式绘出向量a和b定义的离散系统h(k)的时域波形；·impz(b,a,n)绘出向量a和b定义的离散系统在0～n（n必须为整数）离散时间范围内单位脉冲响应h(k)的时域波形；·impz(b,a,n1:n2)绘出向量a和b定义的离散系统在n1～n2（n1、n2必须为整数，且n1<n2）离散时间范围内单位脉冲响应h(k)的时域波形；·y=impz(b,a,n1:n2)不绘出系统的h(k)的时域波形，而是求出向量a和b定义的离散系统在n1～n2（n1、n2必须为整数，且n1<n2）离散时间范围内单位脉冲响应h(k)的数值解。2、stepz函数功能：求离散系统的单位阶跃响应格式：stepz(a,b,c,d)、stepz(num,den)3、dinitial函数功能：求离散系统的零输入响应格式：dinitial(a,b,c,d,x0)、dinitial(num,den,x0)4、dlsim函数功能：求任意输入的离散系统进行仿真格式：dlsim(a,b,c,d,u)、dlsim(num,den,u)5、fliter函数功能：求任意输入的离散系统进行仿真格式：filter(a,b,c,d,u)、filter(num,den,u)设系统差分方程为,求该系统对信号的响应 B=1;A=[1,-0.8]; N=0:31;x=0.8.^n; y=filter(B,A,x); subplot(2,1,1);stem(x) subplot(2,1,2);stem(y)4.3信号与系统的频域分析 4.3.1信号的傅里叶分析 4.3.2连续时间系统的S域分析 4.3.3离散时间系统的Z域分析4.3.1信号的傅里叶分析Matlab的符号运算工具箱（SymbolicMathToolbox）提供了能直接求解傅里叶变换和逆变换的符号运算函数fourier()和ifourier()。两函数的调用格式如下。（1）傅里叶变换在Matlab中，傅里变换变换由函数fourier()实现。fourier()有三种调用格式：①F=fourier(f)求时间函数f(t)的傅里叶变换，返回函数F的自变量默认为w，②F=fourier(f,v)求时间函数f(t)的傅里叶变换，返回函数F的自变量为v，③F=fourier(f,u,v)对自变量为u的函数f(u)求傅里叶变换，返回函数F的自变量为v， （2）傅里叶逆变换 在Matlab中，傅里变换逆变换由函数ifourier()实现。与函数fourier()相类似，ifourier()也有三种调用格式： ①f=ifourier(F) ②f=ifourier(F,u) ③f=ifourier(F,v,u)求单边指数函数的傅里叶变换，画出其幅频特性和相频特性图 symstwf f=exp(-2*t)*sym('Heaviside(t)'); F=fourier(f) subplot(3,1,1);ezplot(f,[0:2,0:1.2]); subplot(3,1,2);ezplot(abs(F),[-10:10]); subplot(3,1,3);ezplot(angle(F),[-10:10]) 严格来说，用数值计算的方法计算连续时间信号的傅里叶变换需要有个限定条件，即信号是时限信号（Timelimitedsignal），也就是当时间|t|大于某个给定时间时其值衰减为零或接近于零，这个条件与前面提到的为什么不能用无限多个谐波分量来合成周期信号的道理是一样的。计算机只能处理有限大小和有限数量的数。 上式用MATLAB表示为： X=x*exp(-j*t’*w)*T采用数值计算算法的理论依据是：若信号为时限信号，当时间间隔T取得足够小时，上式可演变为：二周期连续信号的频域分析-傅里叶级数往往只需对其在一个周期内进行分析即可假设需要计算的谐波次数为N 给定一个周期为T1=2s的连续时间周期方波信号，其一个周期内的数学表达式为 试求其傅里叶级数。3时间离散、连续频率－序列傅里叶变换4离散时间、离散频率－离散傅里叶级数5离散时间、离散频率－离散傅里叶变换(DFT)正变换：逆变换：正变换：逆变换：正变换：逆变换：DFS、DFTfunction[n,X]=dft(x,N)n=[0:N-1];k=[0:N-1];WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;x=[xzeros(1,N-length(x))];X=x*WNnk;1．一维快速正傅里叶变换函数fft格式：X=fft(x,N)功能：采用FFT算法计算序列向量x的N点DFT变换，当N缺省时，fft函数自动按x的长度计算DFT。当N为2整数次幂时，fft按基-2算法计算，否则用混合算法。2．一维快速逆傅里叶变换函数ifft格式：x=ifft(X,N)功能：采用FFT算法计算序列向量X的N点IDFT变换。[例4-36]用快速傅里叶变换FFT计算下面两个序列的卷积。,并测试直接卷积和快速卷积的时间。图4.35快速卷积框图六．快速傅里叶变换MATLAB程序(部分)：%线性卷积xn=sin(0.4*[1:15]);　%对序列x(n)赋值,M=15hn=0.9.^(1:20);　　%对序列h(n)赋值,N=20yn=conv(xn,hn); %直接调用函数conv计算卷积%园周卷积L=pow2(nextpow2(M+N-1));Xk=fft(xn,L); Hk=fft(hn,L); Yk=Xk.*Hk;yn=ifft(Yk,L);图4.36x(n),h(n)及其线性卷积波形4.3.2连续时间系统的S域分析 Matlab的符号数学工具箱（SymbolicMathToolbox）提供了能直接求解拉普拉斯变换和逆变换的符号运算函数laplace()和ilaplace()。1.连续时间信号的拉普拉斯变换 求以下函数的拉普拉斯变换。symstf1=sym('exp(-2*t)*Heaviside(t)');F1=laplace(f1)%求f1(t)的拉普拉斯变换f2=sym('t*exp(-t)*Heaviside(t)');F2=laplace(f2)2.连续系统的复频域分析 系统的频率响应函数是一个非常重要的参数，Matlab工具箱中提供了专门用于分析连续系统频率响应的freqs()函数。利用freqs()函数可求出系统频率响应的数值解，也可绘出其幅频特性和相频特性曲线。 freqs()函数的调用形式有如下四种。 （1）H＝freqs(b,a,w1:dw:w2) 可求得指定频率范围（w1~w2）内相应频点处系统频率响应的样值。其中，w1、w2分别为频率起始值和终止值，dw为频率取样间隔。 （2）[H,w]＝freqs(b,a) 将计算默认频率范围内200个频点上系统频率响应的样值，并赋值给返回向量H，200个频点则记录在向量w中。 （3）[H,w]＝freqs(b,a,n) 将计算默认频率范围内n个频点上系统频率响应的样值，并赋值给返回向量H，n个频点则记录在向量w中。 （4）freqs(b,a) 该调用方式将绘出系统的幅频特性和相频特性曲线。已知一因果系统的系统函数为。试分析其频率特性。 a=[12321]; b=[10-4]; p=roots(a) pxm=max(real(p)); ifpxm>=0 '系统不稳定' else freqs(b,a) end4.3.3离散时间系统的Z域分析 1.离散时间信号的Z变换 MATLAB的符号数学工具箱(SymbolicMathTools)提供了计算z正变换的函数ztrans和计算逆z变换的函数iztrans。其调用形式为： F=ztrans(f)%求符号函数f的z变换，返回函数的自变量为z； F=ztrans(f,w)%求符号函数f的z变换，返回函数的自变量为w； F=ztrans(f,k,w)%对自变量为k的符号函数f求z变换，返回函数的自变量为w。 f=iztrans(F)%对自变量为z的符号函数F求逆z变换，返回函数的自变量为n； f=iztrans(F,k)%对自变量为z的符号函数F求逆z变换，返回函数的自变量为k； f=iztrans(F,w,k)%对自变量为w的符号函数F求逆z变换，返回函数的自变量为k。 2、离散时间系统的Z域分析 在MATLAB中，利用freqz()函数可方便地求得系统的频率响应。调用格式如下： freqz(b，a) 该调用方式将绘制系统在0~p范围内的幅频特性和相频特性图，其中，b、a分别为系统函数分子、分母多项式的系数向量。 freqz(b，a，’whole’) 该调用方式将绘制系统在0~2p范围内的幅频特性和相频特性图。 freqz(b，a，N) 该调用方式将绘制系统在0~p范围内N个频率等分点的幅频特性和相频特性图，N的缺省值为512； freqz(b，a，N，’whole’) 该调用方式将绘制系统在0~2p范围内N个频率等分点的幅频特性和相频特性图。a.画出零极点图b.求系统响应c.求系统的幅频特性和相频特性 b=[0121]; a=[1-0.5-0.0050.3]; subplot(2,2,1),zplane(b,a);%绘制系统的零极点图 title('系统的零极点图'); [hnn]=impz(b,a,16);%用长除法求逆z变换h(n) subplot(2,2,2),stem(n,hn,'filled');%绘制单位响应h(n)的波形 title('单位响应h(n)');gridon;axistight; [hw]=freqz(b,a,16);%计算频率响应 subplot(2,2,3),plot(w,abs(h)');%绘制幅频特性曲线 title('幅频特性曲线');gridon; subplot(2,2,4),plot(w,angle(h));%绘制相频特性曲线 title('相频特性曲线');gridon;*