数字信号处理实验四.doc

实 验 报 告

实验名称_____离散系统分析______

课程名称____数字信号处理________

院 系 部:电气与电子工程 专业班级：信息1002
学生姓名：王萌 学 号: 11012000219
同 组 人: 实验台号:
指导教师：范杰清 成 绩：
实验日期:

华北电力大学
实验四 离散系统分析
一、实验目的
深刻理解离散时间系统的系统函数在分析离散系统的时域特性、频域特性以及稳定性中的重要作用及意义，熟练掌握利用MATLAB分析离散系统的时域响应、频响特性和零极点的方法。掌握利用DTFT和DFT确定系统特性的原理和方法。
二、 实验原理
MATLAB提供了许多可用于分析线性时不变连续系统的函数，主要包含有系统函数、系统时域响应、系统频域响应等分析函数。
1. 离散系统的时域响应
在调用MATLAB函数时，需要利用描述该离散系统的系数函数。对差分方程进行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)：
计算系统的n点频率响应H，w为频率点向量。
H=freqz(b, a, w) ：计算系统在指定频率点向量w上的频响；
freqz(b,a)： 绘制频率响应曲线。
其中：b和a分别为系统函数H(z)的分子分母系数矩阵；
4．利用DTFT和DFT确定离散系统的特性
在很多情况下，需要根据LTI系统的输入和输出对系统进行辨识，即通过测量系统在已知输入x[k]激励下的响应y[k]来确定系统的特性。
若系统的脉冲响应为h[k]，由于存在y[k]=x[k]* h[k]，因而可在时域通过解卷积方法求解h[k]。
在实际应用中，进行信号解卷积比较困难。因此，通常从频域来分析系统，这样就可以将时域的卷积转变为频域的乘积，从而通过分析系统输入序列和输出序列的频率特性来确定系统的频率特性 ，再由 得到系统的脉冲响应h[k]。
若该LTI系统输入x[k]的DTFT为 ,系统输出y[k]的DTFT为 ，则系统的频率特性可表示为：

有限长序列的DTFT可以利用FFT计算出其在区间 内的N个等间隔频率点上的样点值：
X = fft(x,N); Y = fft(y,N);
再利用 H = Y./X 和 h = ifft(H,N)
可以得到系统的单位脉冲响应h[k]。
三、 实验内容
1.已知某LTI系统的差分方程为：

（1）初始状态 ，输入 计算系统的完全响应。
（2）当以下三个信号分别通过系统时，分别计算离散系统的
零状态响应：
（3）该系统具有什么特性？
>> 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);
>> 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)

(1)

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.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.0033 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 through 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 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,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.5616

2.已知某因果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('h[k]')

(2)>>a=[1,-1.035,0.8264,-0.26050.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

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

Columns 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.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.0002 1.0002 1.0002 1.0002

3.已知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)
subplot(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，-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.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,-6.3322,4.9416,-2.648,-0.7564,5.4872,-11.7557,19.7533,-29.6298,41.4666,-55.2433,70.7979,-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,-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命令做出其零极点分布图，利用命令分析该滤波器的幅频特性和相频特性，比较零极点分布与滤波器频率特性的关系。
（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 mtlb
N=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,150);
Y=X.*H;
y=ifft(Y);
subplot(2,3,5);
plot((1:N),y);
title('´¦Àíºó');

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

1.系统函数的零极点对系统频率特性有何影响？
答：系统零点越接近1，则幅频响应的波谷越低，即滤波器的滤波衰减最低点越低；
系统极点越接近1，则幅频响应的波峰越高，即滤波器的滤波通带最高点越高。
极点主要影响频率响应的峰值，极点愈靠近单位圆，峰值愈尖锐；零点主要影响频率特性的谷值，零点愈靠近单位圆，谷值愈深，当零点在单位圆上时，频率特性为零，一个传递函数有几个极点幅度响应就有几个峰值，对应出现一些谷值。
2.对于因果稳定实系数的低通、高通、带通、带阻数字滤波器，极点分布有何特点？
答：因为是因果稳定系统，所以极点都在单位圆内。
3.离散系统的系统函数的零极点对系统脉冲响应有何影响？
答：系统函数的极点位置决定序列包络的变化趋势和变化频率，极点的半径决定了序列包络的变化趋势，而极点的幅角决定序列包络的变化频率，而零点位置只影响冲激响应的幅度大小和相位。
4.若某因果系统不稳定，有哪些主要措施可使之稳定？
答：改变参数，让极点在单位圆内。
5.从频域利用DFT确定离散LTI系统的特性，一般会产生哪些误差，如何改善？
答：频谱混叠，对于带限连续信号，只要提高抽样频率使之满足时域抽样定理；对于非带限信号，更具实际情况对其进行低通滤波，使之成为带限信号；频谱泄露，时域加窗使之成为有限长序列；栅栏现象，在序列后补零，构成新序列后再求频谱
6.若使用DFT对连续LTI系统进行辨识，需要解决哪些问题？
答：频谱混叠的问题，对于带限连续信号，要提高抽样频率使之满足时域抽样定理；对于非带限信号，更具实际情况对其进行低通滤波，使之成为带限信号；