数字信号处理实验二系统响应及系统稳定性

（1）试用Matlab命令求解一下离散时间系统的单位取样响应，并判断系统的稳定性。

    1）3y(n) + 4y(n-1) + y(n-2) = x(n) + x(n-1)

    2)2.5y(n) + 6y(n-1) + 10y(n-2) = x(n)

（2）已知某系统的单位取样响应为，试用 MATLAB求当激励信号为x(n) = u(n) – u(n-5)时，系统的零状态响应。

n = 0:30;
x = (n==0);   %单位序列
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%（1）
figure(1);
subplot(2,2,1);
% 3y1(n)+4y1(n-1)+y1(n-2)=x1(n)+x1(n-1)
a1 = [3 4 1];
b1 = [1 1];
h1 = filter(b1,a1,x);  %使用filter求解单位取样响应
% p = roots(a1);          %求解方程的根，可以根据根的模判断
% impz(b1,a1,30);           % 离散时间系统的单位取样响应
stem(n,h1,'fill');
grid on;
title('1)系统单位取样响应h(n)');

subplot(2,2,2);
zplane(b1,a1);          % 绘制零极点图，可通过极点判断系统是否稳定

subplot(2,2,3);
% 2.5y2(n)+6y2(n-1)+10y(n-2)=x(n)
a2 = [2.5 6 10];
b2 = [1];
% p = roots(a2);          %求解方程的根，可以根据根的模判断
% zplane(b2,a2);          % 绘制零极点图，可通过极点判断系统是否稳定
impz(b2,a2,30);  %使用impz求取单位取样响应
grid on;
title('2)系统单位取样响应h(n)');

subplot(2,2,4);
zplane(b2,a2);          % 绘制零极点图，可通过极点判断系统是否稳定

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%（2）
figure(2);
nh = -1:10;     %将h(n)的范围扩大2个单位
nx = -1:5;      %将x(n)的范围扩大2个单位
hn = (7/8).^nh.*(uDT(nh)-uDT(nh-10));
xn = uDT(nx)-uDT(nx-5);
ny1 = nx(1)+nh(1);      %卷积的起点
ny2 = nx(end)+nh(end);      %卷积的终点
ny = ny1:ny2;           %卷积结果的范围
y = conv(hn,xn);        %卷积的结果
stem(ny,y,'fill','linewidth',1);
xlabel('n');
ylabel('h(n)*x(n)');
title('（2）系统的零状态响应');
grid on;
axis([-3 16 0 5]);
                  

（3）用数字图像处理中的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_colour.bmp');
 % 将图片转换为灰度图
 B = rgb2gray(B);
 
 subplot(2,2,1);
 imshow(B);
 title('原始图像');
 
 img_sobelx = conv2(Gx,B);
 subplot(2,2,2);
 imshow(img_sobelx);
 title('Gx处理过的图像');
 
img_sobely = conv2(Gy,B);
subplot(2,2,3);
imshow(img_sobely);
title('Gy处理过的图像');

img_sobelxy = conv2(Gy,img_sobelx);
subplot(2,2,4);
imshow(img_sobelxy);
title('先Gx后Gy处理过的图像');

思考题（1）Matlab的工具箱函数conv，能用于计算两个有限长序列之间的卷积，但 conv函数假定这两个序列都从n=0开始。试编写M文件计算x(n)=[3,11,7,0,-1,4,2],-3≤n≤3和h(n)=[2,3,0,-5,2,1], -1≤n≤4之间的卷积，并绘制y(n)的波形图。

figure(1);
x = [3 11 7 0 -1 4 2];
nx = -3:3;          %x(n)的范围从-3到3
h = [2 3 0 -5 2 1];
hn = -1:4;          %h(n)的范围从-1到4
y = conv(x,h);
yn = nx(1)+hn(1):nx(end)+hn(end);%卷积的起点是-4，卷积的终点是7
stem(yn,y,'fill');
grid on;
axis([-5 8 -60 60])

思考题（2）在数字图像处理边缘检测技术中，除了 Sobel算子外，还有哪些边缘检测算子？通过查找资料，选择某种边缘检测算子，在Matlab平台上编程实现对lena或其他图像进行边缘检测，显示出原图和边缘检测后的图像。

 % 读取图像文件
 B = imread('lena_colour.bmp');
 % %将图片转换为灰度图
 B = rgb2gray(B);
 
 subplot(2,2,1);
 imshow(B);
 title('原始图像');
 
 img_sobel_both = edge(B,'prewitt',[],'both');  %prewitt算子进行边缘检测，方向为水平和垂直
 subplot(2,2,2);
 imshow(img_sobel_both);
 title('prewitt算子水平和垂直方向处理过的图像');
 
img_sobel_horizontal = edge(B,'prewitt',[],'horizontal');    %prewitt算子进行边缘检测，方向为水平
subplot(2,2,3);
imshow(img_sobel_horizontal);
title('prewitt算子水平方向处理过的图像');

img_sobel_vertical = edge(B,'prewitt',[],'vertical');    %prewitt算子进行边缘检测，方向为垂直
subplot(2,2,4);
imshow(img_sobel_vertical);
title('prewitt算子垂直方向处理过的图像');