实验五 离散时间系统时域分析的MATLAB实现2

实验十一离散信号时域分析的MATLAB实现―、实验目的熟悉MATLAB编程方法、常用语句和可视化绘图技术;掌握序列时域运算的MATLAB编程方法。掌握离散信号时域分析的MATLAB方法。二、实验原理1.离散时间系统LTI离散系统中，其输入和输出的关系由差分方程描述:^^ay(k+i)=Ebf(k+j)(前向差分方程)TOC\o"1-5"\h\zZji=0j=0^^ay(k-i)=£bf(k—n+j)(后向差分方程)Zji=0j=0当系统的输入为单位序列邮)时产生的零状态响应称为系统的单位函数响应，用h(k)表示。当输入为双k)时产生的零状态响应称为系统的单位阶跃应，记为：g(k)，如下图所示。g曲f奇)——或”kLTI离散系统kLTI离散系统►如果系统输入为e(k)，冲激响应为h(k)，系统的零状态响应为y(k)，则有：y(k)=h(k)*f(k)。与连续系统的单位冲激响应h(t)相类似，离散系统的单位函数响应h(k)也包含了系统的固有特性，与输入序列无关。我们只要知道了系统的单位函数响应，即可求得系统在不同激励信号作用下产生的响应。因此，求解系统的单位函数响应h(k)对我们进行离散系统的分析也同样具有非常重要的意义。MATLAB中为用户提供了专门用于求解离散系统单位函数响应，并绘制其时域波形的函数impz()。同样也提供了求离散系统响应的专用函数filter()，该函数能求出由差分方程所描述的离散系统在指定时间范围内的输入序列作用时，产生的响应序列的数值解。当系统初值不为零时，可以使用dlsim()函数求出离散系统的全响应，其调用方法与前面连续系统的lsim()函数相似。另外，求解离散系统阶跃响应可以通过如下两种方法实现：一种是直接调用专用函数dstep(),其调用方法与求解连续系统阶跃响应的专用函数step()的调用方法相似；另一种方法是利用求解离散系统零状态响应的专用函数filter。，只要将其中的激励信号看成是单位阶跃信号s(k)即可。函数的调用格式分别如下：⑴impz()函数impz(b,a)以默认方式绘出由向量a和b所定义的离散系统单位函数响应的时域波形。impz(b,a,n)绘出由向量a和b所定义的离散系统在0〜n(n必须为整数)的离散时间范围内单位函数响应的时域波形。impz(b,a,n1:n2)绘出由向量a和b所定义的离散系统在n1〜n2(n1、n2必须为整数)的离散时间范围内单位函数响应的时域波形。y=impz(b,a,n1:n2)求出由向量a和b所定义的离散系统在n1~n2(n1、n2必须为整数)的离散时间范围内单位函数响应的数值解，但不绘出波形。⑵filter()函数filter(b,a,x)其中a和b与前面相同，x是包含输入序列非零样值点的的行向量。此命令将求出系统在与x的取样时间点相同的输出序列样值。例题：已知描述离散系统的差分方程为：y(k)-0.25y(k-1)+0.5y(k-2)=f(k)+f(k-1)，且已知系统输入序列为f(k)=(1)ke(k)，2求出系统的单位函数响应所k)在-3~10离散时间范围内响应波形。求出系统零状态响应在0~15区间上的样值；并画出输入序列的时域波形以及系统零状态响应的波形分析：①求系统的单位函数响应的MATLAB程序：a=[1,-0.25,0.5];b=[1,1,0];impz(b,a,-3:10),title(单位响应')％绘出单位函数响应在-3~10区间上的波形运行结果如图1。②求零状态响应的MATLAB程序：a=[1,-0.25,0.5];b=[1,1,0]k=0:15;%定义输入序列取值范围x=(1/2).Ak;%定义输入序列表达式y=filter(b,a,x)%求解零状态响应样值subplot(2,1,1),stem(k,x)%绘制输入序列的波形title(输入序列')subplot(2,1,2),stem(k,y)%绘制零状态响应的波形title(输出序列')运行结果如下：y=Columns1through101.00001.75000.6875-0.3281-0.23830.19820.2156-0.0218-0.1015-0.0086Columns11through160.05150.0187-0.0204-0.01410.00690.0088图1.①运行结果图2.②运行结果2.序列的卷积运算序列f1(k)、人(k)的卷积和运算f(k)=f1(k)*/2(k)，可由MATLAB的conv()函数实现，调用格式为：f=conv(f1,f2)[1,-2<k<20,0<k<3如：已知序列：俨)="其他，匕(幻=|0,其他,运行如下M文件可求其卷积和：k1=-2:2;f1=ones(1,length(k1));k2=0:3;f2=2.Ak2;f=conv(f1,f2)结果为：f=137151514128可见，conv()函数不需要给定f](k)、人(k)的非零样值的时间序号，也不返回卷积和序列f(k)的时间序号；此外，conv()假定f](k)、人(k)都是从k=0开始，这就限制了它的应用范围。因此，要对从任意k值开始的序列进行卷积和运算，同时正确标识出函数conv()的计算结果各量f，还须构造序列f](k)、人(k)和犬心的对应序号向量。下面是求序列f](k)、f2(k)卷积和的实用函数dconv()，它可实现序号向量的返回。function[f,k]=dconv(f1,k1,f2,k2)%求卷积和:f(k)=f1(k)*f2(k)f=conv(f1,f2)k0=k1(1)+k2(1);%计算序列f非零样值的起点位置k0k3=length(k1)+length(k2)-2;%计算序列f非零样值的宽度k=k0:k0+k3;%确定序列f非零样值的序号向量

subplot(2,2,1);stem(k1,f1,'fill');title('f1(k)');xlabel('k');subplot(2,2,2);stem(k2,f2,'fill');title('f2(k)');xlabel('k');subplot(2,2,3);stem(k,f,'fill');title('f(k)=f1(k)*f2(k)');xlabel('k');h=get(gca,'position');h(3)=2.5*h(3);set(gca,'position',h)例3.已知序列：f(k)=0.5k,-例3.已知序列：f(k)=0.5k,-2<k<20,其他0,其他图3-2.00004.0000解：运行如下M文件k1=-2:2;f1=0.5.*k1;%定义序列f1(k)k2=0:3;f2=2.Ak2;%定义序列f2(k)[f,k]=dconv(f1,k1,f2,k2);%图3-2.00004.0000图形结果则如图3所示，文本结果如下：f=Columns1through6-1.0000-2.5000-5.0000-9.5000Columns7through88.00008.00003.离散信号相关函数的计算在数字信号处理中，广泛用到信号的相关运算。两个序列的互相关用来度量这两个序列的相似程度。给定两个长度相同、能量受限的序列x(k)和y(k)，它们的互相关函数是另一个序列，定义为：R(n)=芝x(k)y(k-n)=x(n)*y(-n)k=-3其中，n称为滞后参数。当x(k)=y(k)时，称为x(k)的自相关函数，它提供了序列位置不同基准情况下自相似程度的度量。由于信号的自相关函数与信号的功率谱密度是一对傅里叶变换对，这为信号的功率谱计算提供了另一重要途径。由互相关的定义式可知，x(k)和y(k)的互相关函数可以由序列x(n)和y(-n)的卷积和求得，实现互相关函数correlation(x,k1,y,k2)运算的程序如下：function[Rxy,n]=correlation(x,k1,y,k2)[y0,k0]=sigfold(y,k2);%y(k)反转，生成y(-k)[Rxy,n]=dconv(x,k1,y0,k0);holdoffsubplot(3,1,1);stem(k1,x,'fill');xlabel('k');ylabel('x(k)');

subplot(3,1,2);stem(k2,y,'fill');xlabel('k');ylabel('y(k)');subplot(3,1,3);stem(n,Rxy,'fill');xlabel('n');ylabel('Rxy(n)');例4.设x(k)=[2,3,5上一1,3,1]为原型序列，y(k)为x(k)加入噪声干扰并移位后的序列:y(k)=x(k-1)+w(k)其中，w(k)为具有零均值和单位方差的高斯序列。计算序列y(k)和x(k)的互相关。解：运行如下M文件，可得如图4的结果。k1=-3:3;x=[2351-131];%生成x(k)[y2,k2]=sigshift(x,k1,1);%生成x(k-1)w=randn(1,length(y2));%模拟噪声信号w(k)kw=k2;[y,k]=sigadd(y2,k2,w,kw);%生成y(k