实验六 离散系统分析一、实验目的深刻理解离散时间系统的系统函数在分析离散系统的时域特性、频域特性以及稳定性中的重要作用及意义，熟练掌握利用MATLAB分析离散系统的时域响应、频响特性和零极点的方法。掌握利用DTFT和DFT确定系统特性的原理和方法。二、实验内容1.已知某LTI系统的差分方程为：（1）初始状态 ，输入 计算系统的完全响应。（2）当以下三个信号分别通过系统时，分别计算离散系统的零状态响应：（3）该系统具有什么特性？解（1）代码如下：a=1,-1.143,0.412;b=0.0675,0.1349,0.0675;N=50;x=ones(1,N);zi=filtic(b,a,1,2);y=filter(b,a,x,zi)stem(y);xlabel(k);title(hk)运行结果：（2）代码如下：a=1,-1.143,0.412;b=0.0675,0.1349,0.0675;N=25;k=1:N;x1=cos(pi/10*k);y1=filter(b,a,x1)stem(y1);运行结果：x2=cos(pi/5*k);y2=filter(b,a,x2)stem(y2);x3=cos(7*pi/10*k);y3=filter(b,a,x3)stem(y3);（3）确定该微分方程的零极点:a=1,-1.143,0.4128; b=0.0675,0.1349,0.0675; z=roots(b) p=roots(a) zplane(b,a)结果如下：由上图可知，系统的零极点都在单位圆内，由此可得，该系统具有稳定性和因果性。2已知某因果LTI离散系统的系统函数为：(1) 计算系统的单位脉冲响应。(2) 当信号通过系统时，计算系统的零状态响应。解：（1）代码如下: N=40; a=1,-1.035,0.8264,-0.2605,0.04033; b=0.03571,0.1428,0.2143,0.1428,0.03571; y=impz(b,a,N); stem(y) xlabel(k); title(hk);运行结果：（2）代码如下：k=0:49; a=1,-1.035,0.8264,-0.2605,0.04033; b=0.03571,0.1428,0.2143,0.1428,0.03571; x=1+cos(pi/4*k)+cos(pi/2*k); y=filter(b,a,x) stem(k,y, fill);grid on; xlabel(k);title(系统零状态响应y(k);运行结果：3已知LTI系统的输入输出序列分别为(a) (b) (1) 利用系统辨识原理从频域分别求解系统的单位脉冲响应.(2) 利用解析方法确定，并求出系统的单位脉冲响应。比较解析方法与系统辨识方法得到的系统单位脉冲响应，分析误差原因。（1）代码如下：N=10;K=0:N-1;aa=1,-0.25;ab=5/4,-7/16;ba=1;bb=1,-1;ay=impz(ab,aa,N);subplot(2,1,1);stem(ay);by=impz(bb,ba,N);subplot(2,1,2);stem(by);运行结果：（2）k=1:100;x1=(1/2).k;y1=(1/4).*(1/2).k+(1/4).k;X1=fft(x1);Y1=fft(y1);H1=Y1/X1;H1=ifft(H1)x2=(1/4).k;y2=(1/4).k-(1/4).(k-1);X2=fft(x2);Y2=fft(y2);H2=Y2/X2;h2=ifft(H2)结果：H1 = 0.6786 + 0.0000ih2 = -3.00004利用load mtlb命令读入一段语音信号得到序列，然后在该段语音信号中加入500Hz的正弦型干扰信号得到信号，利用FFT分析其频谱。(1) 下列数字滤波器能够滤除信号中500Hz正弦型干扰信号利用zplane命令做出其零极点分布图，利用freqz命令分析该滤波器的幅频特性和相频特性，比较零极点分布与滤波器频率特性的关系。(2) 利用该数字滤波器滤除信号中的噪声，利用FFT观察其频谱，利用sound函数播放处理前后的信号，比较处理前后的效果。解：（1） 代码如下：a=1,-3.594,5.17,-3.494,0.945;b=0.6877,-2.509,3.664,-2.509,0.6877;z=roots(b);p=roots(a);subplot(2,2,1);zplane(b,a);H,w=freqz(b,a);subplot(2,2,2);plot(w,abs(H);xlabel(w(rad);ylabel(H0);title(幅频特性);subplot(2,2,3);plot(w,angle(H);xlabel(w(rad);ylabel( );title(相频特性);运行结果：（2）代码如下： load mtlbN=150;x=mtlb(1:N);k=1:N;subplot(1,2,1);plot(k,x);title(处理前);X=fft(x,150);h=ifft(H);H=fft(H,150);Y=X.*H;y=ifft(Y);subplot(1,2,2);plot(1:N),y);title(处理后);运行结果：