实验二：系统响应及系统稳定性

实验内容

（1）离散时间系统的单位取样响应

试用 Matlab 命令求解以下离散时间系统的单位取样响应，并判断系统的稳定性。
1）
2）
figure;
% （1）
subplot(2, 1, 1);
a = [3, 4, 1]; % 差分方程右端的系数向量
b = [1, 1]; % 差分方程左端的系数向量
n = 0: 10;
% x = (0.5) .^ n; % 激励信号
x = (n == 0); % 冲激
y = filter(b, a, x); % 零状态响应
stem(n, y, 'fill');
title({ ...
'单位取样响应 y(n)'; ...
'3y(n)+ 4y(n-1) + y(n-2) = x(n)+ x(n-1)' ...
});
xlabel('n');
ylabel('y');
grid on;
% （2）
subplot(2, 1, 2);
a = [5/2, 6, 10]; % 差分方程右端的系数向量
b = [1]; % 差分方程左端的系数向量
n = 0: 10;
% x = (0.5) .^ n; % 激励信号
x = (n == 0); % 冲激
y = filter(b, a, x); % 零状态响应
stem(n, y, 'fill');
title({ ...
'单位取样响应 y(n)'; ...
'5/2y(n)+ 6y(n-1) + 10y(n-2) = x(n)' ...
});
xlabel('n');
ylabel('y');
grid on;

（2）离散时间系统的零状态响应

单位取样响应：
求系统的激励信号为：的零状态响应
nx = -3: 11;
x = uDT(nx) - uDT(nx - 5); % 激励信号
nh = -5: 11;
h = (7/8) .^ nh .* (uDT(nh) - uDT(nh - 10)); % 单位取样响应
ny1 = nx(1) + nh(1);
ny2 = nx(end) + nh(end);
ny = ny1: ny2;
y = conv(x, h); % 零状态响应
limits = [-10, 24, 0, 5];
subplot(311)
stem(nx, x, 'fill');
grid on;
xlabel('n');
title('激励信号 x(n)')
axis(limits);
subplot(312);
stem(nh, h, 'fill');
grid on;
xlabel('n');
title('系统的单位取样响应 h(n)');
axis(limits);
subplot(313);
stem(ny, y, 'fill');
grid on;
xlabel('n');
title('系统的零状态响应 y(n) = x(n)*h(n)');
axis(limits);

（3）用数字图像处理中的 Sobel 算子，调用 conv2 函数对 lena 图像进行滤波。

用数字图像处理中的 Sobel 算子，调用 conv2 函数对 lena 图像进行滤波。设两个 3×3 的 Sobe 矩阵算子 Gx=[-1 0 1;-2 0 2;-1 0 1]，Gy=[1 2 1;0 0 0;-1 -2 -1]。
1）读取 lena 图像数据，保存到 B 矩阵中；
2）分如下三种情况对 B 进行卷积滤波
① 采用 Gx 与 B 进行卷积；
② 采用 Gy 与的 B 进行卷积；
③ 先采用 Gx 与 B 进行卷积，然后再采用 Gy 与 B 进行卷积；
④ 在一个 figure 中，分四个子图显示出：原始图像，经过 Gx 卷积滤波后图像、经过 Gy 卷积滤波后图像，先后采用 Gx 和 Gy 与的 B 进行卷积滤波后的图像，观察滤波后的图像的效果。
%% 灰度图卷积
Gx = [-1 0 1; -2 0 2; -1 0 1];
Gy = [1 2 1; 0 0 0; -1 -2 -1];
B = imread('lena.bmp'); % 读取图像
res1 = conv2(B, Gx);
res2 = conv2(B, Gy);
res3 = res1 + res2; % 这个效果比下面那个好
% res3 = conv2(res1, Gy); % 这个不行
figure;
subplot(2, 2, 1);
imshow(B);
title('原始图像');
subplot(2, 2, 2);
imshow(res1);
title('Gx 卷积滤波后图像');
subplot(2, 2, 3);
imshow(res2);
title('Gy 卷积滤波后图像');
subplot(2, 2, 4);
imshow(res3);
title('Gx 和 Gy 卷积滤波后图像');
%% 彩图卷积
Gx = [-1 0 1; -2 0 2; -1 0 1];
Gy = [1 2 1; 0 0 0; -1 -2 -1];
B = imread('lena_colour.bmp'); % 读取图像
[row, col, chan] = size(B);
figure;
subplot(2, 2, 1);
imshow(B);
title('原始图像');
subplot(2, 2, 2);
fig1 = zeros(row, col, 3);
for i = 1:3
fig1(:,:,i) = conv2(B(:,:,i), Gx, 'same');
end
imshow(fig1);
title('Gx 卷积滤波后图像');
subplot(2, 2, 3);
fig2 = zeros(row, col, 3);
for i = 1:3
fig2(:,:,i) = conv2(B(:,:,i), Gy, 'same');
end
imshow(fig2);
title('Gy 卷积滤波后图像');
subplot(2, 2, 4);
fig3 = fig1 + fig2;
imshow(fig3);
title('Gx 和 Gy 卷积滤波后图像');

思考题

计算离散时间系统的卷积和、绘制卷积波形图

（1） Matlab 的工具箱函数 conv，能用于计算两个有限长序列之间的卷积，但 conv 函 数 假 定 这 两 个 序 列 都 从 n=0 开 始 。 试 编 写 M 文 件 计 算 x(n) 和 h(n) 之 间 的 卷积，并绘制 y(n) 的波形图。
（2） 在数字图像处理边缘检测技术中，除了 Sobel 算子外，还有哪些边缘检测算子？通过查找资料，选择某种边缘检测算子，在 Matlab 平台上编程实现对 lena 或其他图像进行边缘检测，显示出原图和边缘检测后的图像。
说明：
（1）计算、绘制 y(n)，编写的 convu() 在最下面
figure;
nx = -3: 3;
x = [3, 11, 7, 0, 1, 4, 2];
nh = -1: 4;
h = [2, 3, 0, 5, 2, 1];
[y, ny] = convu(h, nh, x, nx);
subplot(3, 1, 1);
stem(nx, x, "filled");
grid on;
xlabel('n');
ylabel('x(n)');
title('x(n) 的波形图');
subplot(3, 1, 2);
stem(nh, h, "filled");
grid on;
xlabel('n');
ylabel('h(n)');
title('h(n) 的波形图');
subplot(3, 1, 3);
stem(ny, y, 'fill');
grid on;
xlabel('n');
ylabel('y(n)');
title('y(n) 的波形图');
（2）显示出原图和边缘检测后的图像
B = imread('lena.bmp'); % 读取图像
figure;
subplot(2, 2, 1);
imshow(B);
title('原始图像');
subplot(2, 2, 2);
fig1 = edge(B, 'canny');
imshow(fig1);
title('canny 算子边缘检测后的图像');
subplot(2, 2, 3);
fig2 = edge(B, 'sobel');
imshow(fig2);
title('sobel 算子边缘检测后的图像');
subplot(2, 2, 4);
fig3 = edge(B, 'log');
imshow(fig3);
title('log 算子边缘检测后的图像');
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