实验三：z 变换及离散 LTI 系统的 z 域分析

实验内容

1、有理函数 z 变换的部分分式展开和

topic = "部分分式展开和";
figure;
B = [2, 16, 44, 56, 32]; % 分子多项式的系数向量
A = [3, 3, -15, 18, -12]; % 分母多项式的系数向量
[R, P, K] = residuez(B, A); % R 为部分分式的系数向量；P 为极点向量；K 为多项式的系数
% topic, R, P, K
disp(topic);
部分分式展开和
disp(R);
-0.0177 + 0.0000i 9.4914 + 0.0000i -3.0702 + 2.3398i -3.0702 - 2.3398i
disp(P);
-3.2361 + 0.0000i 1.2361 + 0.0000i 0.5000 + 0.8660i 0.5000 - 0.8660i
disp(K);
-2.6667

2、画出系统函数零极点分布图，判断系统的稳定性

（1）
（2）
注意：
tf2zp() 函数输入的 B 和 A 长度要一致，不足的在前面补 0，即 z 的高次系数为 0
topic = "零极点";
%% （1）
h = figure;
% subplot(1, 2, 1);
B = [0, 2, -1.6, -0.9];
A = [1, -2.5, 1.96, -0.48];
[R, P, K] = tf2zp(B, A); % R 为零点；P 为极点；K 为增益
% topic, R, P, K
disp(topic);
零极点
disp(R);
1.1810 -0.3810
disp(P);
1.2000 0.8000 0.5000
disp(K);
2
zplane(B, A);
grid on;
legend('零点', '极点');
title('零极点分布图');
saveas(h, '问题2-零极点分布图1', 'svg');
% 存在极点在 z 域单位圆外，故不稳定
%% （2）
h = figure;
% subplot(1, 2, 2);
B = [0, 0, 0, 1, -1];
A = [1, -0.9, -0.65, 0.873, 0];
[R, P, K] = tf2zp(B, A); % R 为零点；P 为极点；K 为增益
% topic, R, P, K
disp(topic);
零极点
disp(R);
1
disp(P);
0.0000 + 0.0000i -0.9000 + 0.0000i 0.9000 + 0.4000i 0.9000 - 0.4000i
disp(K);
1
zplane(B, A);
grid on;
legend('零点', '极点');
title('零极点分布图');
saveas(h, '问题2-零极点分布图2', 'svg');
% 所有极点在 z 域单位圆内，故稳定
% 稳定：系统的极点全部在单位圆内
% 不稳定：系统的极点有在单位圆外

3、绘制系统的频率响应曲线

h = figure;
B = [1, 0, 0];
A = [1, -3/4, 1/8];
[H, W] = freqz(B, A, 400, "whole");
Hm = abs(H);
Hp = angle(H);
subplot(2, 1, 1);
plot(W / pi, Hm);
grid on;
xlabel('\omega(\pi rad/s)');
ylabel('Magnitude');
title('离散系统的幅频特性曲线');
subplot(2, 1, 2);
plot(W / pi, Hp);
grid on;
xlabel('\omega(\pi rad/s)');
ylabel('Phase');
title('离散系统的相频特性曲线');
saveas(h, '问题3-系统的频率响应曲线', 'svg');

4、画零极图、判断稳定性、画幅/相频特性曲线、判断系统

已知系统的差分方程：
（1）画出该系统的零极点分布图，判断系统的稳定性
（2）画出系统在 范围内的幅频特性曲线和相频特性曲线
（3）查找资料说明该系统的功能
答：
（1）
（2）如图所示
（3）功能：梳妆滤波器
%% （1）
topic = "零极点";
h = figure;
% subplot(1, 2, 1);
B = [1, 0, 0, 0, 0, 0, 0, 0, -1];
A = [1, 0, 0, 0, 0, 0, 0, 0, -0.9];
[R, P, K] = tf2zp(B, A); % R 为零点；P 为极点；K 为
% topic, R, P, K;
disp(topic);
零极点
disp(R);
-1.0000 + 0.0000i -0.7071 + 0.7071i -0.7071 - 0.7071i 0.0000 + 1.0000i 0.0000 - 1.0000i 1.0000 + 0.0000i 0.7071 + 0.7071i 0.7071 - 0.7071i
disp(P);
-0.9869 + 0.0000i -0.6979 + 0.6979i -0.6979 - 0.6979i -0.0000 + 0.9869i -0.0000 - 0.9869i 0.9869 + 0.0000i 0.6979 + 0.6979i 0.6979 - 0.6979i
disp(K);
1
zplane(B, A);
grid on;
% axis(limits);
legend('零点', '极点');
title('零极点分布图');
saveas(h, '问题4-零极图', 'svg');
% 稳定
%% （2）
h = figure;
[H, W] = freqz(B, A, 400, "whole");
Hm = abs(H);
Hp = angle(H);
subplot(2, 1, 1);
plot(W / pi, Hm);
grid on;
xlabel('\omega(\pi rad/s)');
ylabel('Magnitude');
title('离散系统的幅频特性曲线');
subplot(2, 1, 2);
plot(W / pi, Hp);
grid on;
xlabel('\omega(\pi rad/s)');
ylabel('Phase');
title('离散系统的相频特性曲线');
saveas(h, '问题4-频率响应曲线', 'svg');
%% （3）功能：梳妆滤波器

思考题

1、分别采用系统
对音频文件motherland.wav进行滤波（可采用实验二的conv函数）
（1）画出滤波前后该音频文的连续时域波形图
（2）分析说明滤波后信号幅度变化的原因
笔记：
高通，声音变尖；低通：声音变低
斜变信号
不稳定
说明：
默认不播放音频，需要可自行注释，部分音频出现“没有波形”属于正常，因为部分音频不全，在全部音频下可以看到波形
[xn, fs] = audioread('motherland.wav'); % [音频数据, 采样频率]
t = (0: length(xn) - 1) / fs;
left = find(t == 1);
right = find(t == 2);
n = right - left + 1;
% sound(xn(left:right)); % 1~2s
%% （1）系统1
B1 = [1, 0];
A1 = [1, 0.8];
hn = filter(B1, A1, xn);
% % 全部音频
% figure;
% subplot(2, 1, 1);
% plot(t, xn);
% grid on;
% xlabel('时间 t');
% title('滤波前音频文件的连续时域波形图');
%
% subplot(2, 1, 2);
% plot(t, hn);
% grid on;
% xlabel('时间 t');
% title('滤波后音频文件的连续时域波形图');
% 部分音频 1-2s
h = figure;
subplot(2, 1, 1);
plot(t(left:right), xn(left:right));
grid on;
xlabel('时间 t');
title('滤波前音频文件的连续时域波形图');
subplot(2, 1, 2);
plot(t(left:right), hn(left:right));
grid on;
xlabel('时间 t');
title('滤波后音频文件的连续时域波形图');
% sound(hn(left:right)); % 播放音频
saveas(h, '思考题1-系统1', 'svg');
%% （2）系统2
B2 = [1, 0];
A2 = [1, -1];
hn = filter(B2, A2, xn);
% % 全部音频
% figure;
% subplot(2, 1, 1);
% plot(t, xn);
% grid on;
% xlabel('时间 t');
% title('滤波前音频文件的连续时域波形图');
%
% subplot(2, 1, 2);
% plot(t, hn);
% grid on;
% xlabel('时间 t');
% title('滤波后音频文件的连续时域波形图');
% 部分音频 1-2s
h = figure;
subplot(2, 1, 1);
plot(t(left:right), xn(left:right));
grid on;
xlabel('时间 t');
title('滤波前音频文件的连续时域波形图');
subplot(2, 1, 2);
plot(t(left:right), hn(left:right));
grid on;
xlabel('时间 t');
title('滤波后音频文件的连续时域波形图');
% sound(hn(left:right)); % 播放音频
saveas(h, '思考题1-系统2', 'svg');
%% （3）系统3
B3 = [1, 0];
A3 = [1, 1.2];
hn = filter(B3, A3, xn);
% % 全部音频
% figure;
% subplot(2, 1, 1);
% plot(t, xn);
% grid on;
% xlabel('时间 t');
% title('滤波前音频文件的连续时域波形图');
%
% subplot(2, 1, 2);
% plot(t, hn);
% grid on;
% xlabel('时间 t');
% title('滤波后音频文件的连续时域波形图');
% 部分音频 1-2s
h = figure;
subplot(2, 1, 1);
plot(t(left:right), xn(left:right));
grid on;
xlabel('时间 t');
title('滤波前音频文件的连续时域波形图');
subplot(2, 1, 2);
plot(t(left:right), hn(left:right));
grid on;
xlabel('时间 t');
xlim([1, 2]);
title('滤波后音频文件的连续时域波形图');
% sound(hn(left:right)); % 播放音频
saveas(h, '思考题1-系统3', 'svg');
2、已知系统的差分方程：
（1）画出该系统的零极点分布图，判断系统的稳定性
（2）画出系统在 范围内的幅频特性曲线和相频特性曲线
（3）查找资料说明该系统的功能
（1）
（2）如图所示
（3）功能：全通滤波器
%% （1）
topic = "零极点";
h = figure;
% subplot(1, 2, 1);
B = [1, -2, 2];
A = [2, -2, 1];
% R 为零点；P 为极点；K 为
[R, P, K] = tf2zp(B, A);
% topic, R, P, K;
disp(topic);
零极点
disp(R);
1.0000 + 1.0000i 1.0000 - 1.0000i
disp(P);
0.5000 + 0.5000i 0.5000 - 0.5000i
disp(K);
0.5000
zplane(B, A);
grid on;
% axis(limits);
legend('零点', '极点');
title('零极点分布图');
% 稳定
saveas(h, '思考题2-零极点分布图', 'svg');
%% （2）
h = figure;
[H, W] = freqz(B, A, 400, "whole");
Hm = abs(H);
Hp = angle(H);
subplot(2, 1, 1);
plot(W / pi, Hm);
grid on;
xlabel('\omega(\pi rad/s)');
ylabel('Magnitude');
title('离散系统的幅频特性曲线');
subplot(2, 1, 2);
plot(W / pi, Hp);
grid on;
xlabel('\omega(\pi rad/s)');
ylabel('Phase');
title('离散系统的相频特性曲线');
saveas(h, '思考题2-频率响应曲线', 'svg');
%% （3）功能：全通滤波器
function [y, ny] = convu(h, nh, x, nx)
ny1 = nh(1) + nx(1);
ny2 = nh(end) + nx(end);
y = conv(h, x);
ny = ny1: ny2;
end