数字信号处理实验四

1、实验报告课程名称：数字信号处理实 验 四：离散系统分析 班 级：通信1403 学生姓名：强亚倩 学 号：1141210319 指导教师：范杰清华北电力大学（北京）一、实验目的 深刻理解离散时间系统的系统函数在分析离散系统的时域特性、频域特性以及稳定性中的重要作用及意义，熟练掌握利用MATLAB分析离散系统的时域响应、频响特性和零极点的方法。掌握利用DTFT和DFT确定系统特性的原理和方法。二、实验原理MATLAB提供了许多可用于分析线性时不变连续系统的函数，主要包含有系统函数、系统时域响应、系统频域响应等分析函数。1. 离散系统的时域响应在调用MATLAB函数时，需要利用描述该离散系统的系数函

2、数。对差分方程进行Z变换即可得系统函数：在MATLAB中可使用向量a和向量b分别保存分母多项式和分子多项式的系数：这些系数均从z0按z的降幂排列。2离散系统的系统函数零极点分析离散LTI系统的系统函数H(z)可以表示为零极点形式：使用MATLAB提供的roots函数计算离散系统的零极点；使用zplane函数绘制离散系统的零极点分布图。注意：在利用这些函数时，要求H(z)的分子多项式和分母多项式的系数的个数相等，若不等则需要补零。3离散系统的频率响应当离散因果LTI系统的系统函数H(z)的极点全部位于z平面单位圆内时，系统的频率响应可由H(z)求出，即 H, w=freqz(b, a, n)：

3、计算系统的n点频率响应H，w为频率点向量。H=freqz(b, a, w) ：计算系统在指定频率点向量w上的频响；freqz(b,a)： 绘制频率响应曲线。其中：b和a分别为系统函数H(z)的分子分母系数矩阵；4利用DTFT和DFT确定离散系统的特性 在很多情况下，需要根据LTI系统的输入和输出对系统进行辨识，即通过测量系统在已知输入xk激励下的响应yk来确定系统的特性。若系统的脉冲响应为hk，由于存在yk=xk* hk，因而可在时域通过解卷积方法求解hk。在实际应用中，进行信号解卷积比较困难。因此，通常从频域来分析系统，这样就可以将时域的卷积转变为频域的乘积，从而通过分析系统输入序列和输出序

4、列的频率特性来确定系统的频率特性 ，再由 得到系统的脉冲响应hk。若该LTI系统输入xk的DTFT为 ,系统输出yk的DTFT为 ，则系统的频率特性可表示为：有限长序列的DTFT可以利用FFT计算出其在区间 内的N个等间隔频率点上的样点值： X = fft(x,N); Y = fft(y,N);再利用 H = Y./X 和 h = ifft(H,N) 可以得到系统的单位脉冲响应hk。三、 实验内容1.已知某LTI系统的差分方程为：（1）初始状态 ，输入 计算系统的完全响应。（2）当以下三个信号分别通过系统时，分别计算离散系统的零状态响应：（3）该系统具有什么特性？(1) N=100;b=0.0

5、675,0.1349,0.0675; a=1,-1.143,0.412;x=ones(1,N);zi=filtic(b,a,1,2);y=filter(b,a,x,zi); N=100;b=0.0675,0.1349,0.0675; a=1,-1.143,0.412;x=ones(1,N);zi=filtic(b,a,1,2);y=filter(b,a,x,zi)y = Columns 1 through 6 0.3865 0.2322 0.3760 0.6041 0.8054 0.9416 Columns 7 through 12 1.0143 1.0413 1.0422 1.0322 1.

6、0202 1.0108 Columns 13 through 18 1.0049 1.0020 1.0012 1.0015 1.0021 1.0027 Columns 19 through 24 1.0031 1.0033 1.0034 1.0035 1.0034 1.0034 Columns 25 through 30 1.0034 1.0034 1.0033 1.0033 1.0033 1.003 Columns 31 through 36 1.0033 1.0033 1.0033 1.0033 1.0033 1.0033 Columns 37 through 42 1.0033 1.00

7、33 1.0033 1.0033 1.0033 1.0033 Columns 43 through 48 1.0033 1.0033 1.0033 1.0033 1.0033 1.0033 Columns 49 through 54 1.0033 1.0033 1.0033 1.0033 1.0033 1.0033 Columns 55 through 60 1.0033 1.0033 1.0033 1.0033 1.0033 1.0033 Columns 61 through 66 1.0033 1.0033 1.0033 1.0033 1.0033 1.0033 Columns 67 th

8、rough 72 1.0033 1.0033 1.0033 1.0033 1.0033 1.0033 Columns 73 through 78 1.0033 1.0033 1.0033 1.0033 1.0033 1.0033 Columns 79 through 84 1.0033 1.0033 1.0033 1.0033 1.0033 1.0033 Columns 85 through 90 1.0033 1.0033 1.0033 1.0033 1.0033 1.0033 Columns 91 through 96 1.0033 1.0033 1.0033 1.0033 1.0033

9、1.0033 Columns 97 through 100 1.0033 1.0033 1.0033 1.0033 b=0.0675,0.1349,0.0675; a=1,-1.143,0.412;k=0,0.001,100;x=cos(pi/10*k);y=filter(b,a,x) y= 0.0675 0.2796 0.5616>> b=0.0675,0.1349,0.0675; a=1,-1.143,0.412;k=0,0.001,100;x=cos(pi/5*k);y=filter(b,a,x)y =0.0675 0.2796 0.5616>> b=0.0675

10、,0.1349,0.0675; a=1,-1.143,0.412;k=0,0.001,100;x=cos(7*pi/10*k);y=filter(b,a,x)y = 0.0675 0.2796 0.56162.已知某因果LTI系统的系统函数为：（1）计算系统的单位脉冲响应。 (2）当信号通过系统时，计算系统的零状态响应。a=1,-1.035,0.8264,-0.2605,0.04033; b=0.03571,0.1428,0.2143,0.1428,0.03571;N=100;y=impz(b,a,N);stem(y)xlabel('k');title('hk'

11、)(2)>>a=1,-1.035,0.8264,-0.2605,0.04033; b=0.03571,0.1428,0.2143,0.1428,0.03571;N=100;u=ones(1,N);x=u+cos(pi/4*N)+cos(pi/2*N);y=filter(b,a,x)结果：y = Columns 1 through 6 0.0357 0.2155 0.5863 0.9737 1.1492 1.1002 Columns 7 through 12 0.9903 0.9472 0.9735 1.0098 1.0187 1.0066 Columns 13 through 18

12、 0.9951 0.9940 0.9989 1.0024 1.0021 1.0002 Columns 19 through 24 0.9993 0.9996 1.0002 1.0005 1.0003 1.0001 Columns 25 through 30 1.0000 1.0001 1.0002 1.0002 1.0002 1.0001 Columns 31 through 36 1.0001 1.0002 1.0002 1.0002 1.0002 1.0002 Columns 37 through 42 1.0002 1.0002 1.0002 1.0002 1.0002 1.0002 C

13、olumns 43 through 48 1.0002 1.0002 1.0002 1.0002 1.0002 1.0002 Columns 49 through 54 1.0002 1.0002 1.0002 1.0002 1.0002 1.0002 Columns 55 through 60 1.0002 1.0002 1.0002 1.0002 1.0002 1.0002 Columns 61 through 66 1.0002 1.0002 1.0002 1.0002 1.0002 1.0002 Columns 67 through 72 1.0002 1.0002 1.0002 1.

14、0002 1.0002 1.0002 Columns 73 through 78 1.0002 1.0002 1.0002 1.0002 1.0002 1.0002 Columns 79 through 84 1.0002 1.0002 1.0002 1.0002 1.0002 1.0002 Columns 85 through 90 1.0002 1.0002 1.0002 1.0002 1.0002 1.0002 Columns 91 through 96 1.0002 1.0002 1.0002 1.0002 1.0002 1.0002 Columns 97 through 100 1.

15、0002 1.0002 1.0002 1.00023.已知LTI系统的输入输出序列分别为（1）利用系统辨识原理从频域分别求解系统的单位脉冲响应 。（2）利用解析方法确定 ，并求出系统的单位脉冲响应 。比较解析方法与系统辨识方法得到的系统单位脉冲响应，分析误差原因。(1)N=10; （2）k=1:100;K=0:N-1; x1=(1/2).k; aa=1,-0.25; y1=(1/4).*(1/2).k+(1/4).k;ab=5/4,-7/16; X1=fft(x1);Y1=fft(y1);ba=1;bb=1,-1; H1=Y1/X1;ay=impz(ab,aa,N); H1=ifft(H1)s

16、ubplot(2,1,1); x2=(1/4).k;stem(ay);by=impz(bb,ba,N); y2=(1/4).k-(1/4).(k-1);subplot(2,1,2); X2=fft(x2);Y2=fft(y2);stem(by); H2=Y2/X2; h2=ifft(H2)4.已知某离散系统的输入输出序列。输入序列：2，0.8333，0.3611，0.162，0.0748，0.0354，0.017，0.0083，0.0041，0.002，0.001，0.0005，0.0002，0.0001，0.0001，后面的数值均趋于0；输出序列：0.0056，-0.0259，0.073，-

17、0.1593，0.297， -0.4974，0.7711，-1.1267，1.5702，-2.1037，2.724，-3.4207，4.174，-4.9528，5.7117，-6.3889，6.9034，-7.1528，7.012，-6.3322，4.9416，-2.648，-0.7564，5.4872，-11.7557，19.7533，-29.6298，41.4666，-55.2433，70.7979，-87.7810（1）绘出输入输出信号的波形。（2）计算该系统的频率响应 ，并绘出其幅频特性。（3）计算该系统的单位脉冲响应 ，并绘出其波形。(1)x=2,0.8333,0.3611,0.16

18、2,0.0748,0.0354,0.017,0.0083,0.0041,0.002,0.001,0.0005,0.0002,0.0001,0.0001;y=0.0056,-0.0259,0.073,-0.1593,0.297,-0.4974,0.7711,-1.1267,1.5702,-2.1037,2.724,-3.4207,4.174,-4.9528,5.7117,-6.3889,6.9034,-7.1528,7.012,-6.3322,4.9416,-2.648,-0.7564,5.4872,-11.7557,19.7533,-29.6298,41.4666,-55.2433,70.79

19、79,-87.7810;subplot(2,1,1);stem(x);subplot(2,1,2);stem(y); (2)x=2,0.8333,0.3611,0.162,0.0748,0.0354,0.017,0.0083,0.0041,0.002,0.001,0.0005,0.0002,0.0001,0.0001;y=0.0056,-0.0259,0.073,-0.1593,0.297,-0.4974,0.7711,-1.1267,1.5702,-2.1037,2.724,-3.4207,4.174,-4.9528,5.7117,-6.3889,6.9034,-7.1528,7.012,-

20、6.3322,4.9416,-2.648,-0.7564,5.4872,-11.7557,19.7533,-29.6298,41.4666,-55.2433,70.7979,-87.7810;X=fft(x,256);Y=fft(y,256);H=Y./X;H0=abs(H);plot(H0);(3)h=ifft(H); plot(h0);5.利用load mtlb命令读入一段语音信号得到序列 ，然后在该段语音信号中加入500Hz的正弦型干扰信号得到信号 ，利用FFT分析其频谱。（1）下列数字滤波器能够滤除信号 中500Hz正弦型干扰信号。 利用zplane命令做出其零极点分布图，利用命令分析

21、该滤波器的幅频特性和相频特性，比较零极点分布与滤波器频率特性的关系。（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'

22、);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(2,3,4);plot(k,x);title('´¦À

23、;íÇ°');X=fft(x,150);h=ifft(H);H=fft(H,150);Y=X.*H;y=ifft(Y);subplot(2,3,5);plot(1:N),y);title('´¦Àíºó');(2) load mtlbN=150;x=mtlb(1:N);k=1:N;subplot(2,3,4);plot(k,x);title('´¦ÀíÇ°');X=fft(x,150);h=ifft(H);H=fft(H,15