数字信号与图像处理实验一：信号处理基础

掌握信号处理的基本思想，理解采样信号的频谱特性,加强信号采样与重建的有关基本概念的理解，深入理解线性时不变系统输出与输入的关系，了解数字信号采样率转换前后信号频谱的特征。

文章目录

数字信号与图像处理实验一：信号处理基础
- 1. 给定序列，绘制其图像并分析其混叠
- - 1.1 用MATLAB实现序列 x ( n ) = 3 c o s ( 0.125 π n + 0.2 π ) + 2 s i n ( 0.25 π n + 0.1 π ) x(n) = 3cos(0.125\pi n+0.2\pi)+2sin(0.25\pi n+0.1\pi) x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π)。
  - 1.2 将1.1中的 x ( n ) x(n) x(n)扩展为以16为周期的函数 x ( n ) 16 = x ( n + 16 ) x(n)_{16}=x(n+16) x(n)16=x(n+16)，绘制出4个周期。
  - 1.3 将1.1中的 x ( n ) x(n) x(n)扩展为以10为周期的函数 x ( n ) 16 = x ( n + 10 ) x(n)_{16}=x(n+10) x(n)16=x(n+10)，绘制出4个周期。
- 2. 设计序列的移位和累加函数
- - 2.1 x 1 ( n ) = 2 x ( n + 2 ) − x ( n − 1 ) − 2 x ( n ) x_1(n)=2x(n+2)-x(n-1)-2x(n) x1(n)=2x(n+2)−x(n−1)−2x(n)
  - 2.2 x 2 ( n ) = n ∑ k = 1 5 x ( n − k ) x_2(n)=n\sum_{k=1}^5x(n-k) x2(n)=n∑k=15x(n−k)
- 3. 绘制因果系统的幅频响应和相频响应
- 4. 验证DTFT的折叠特性
- 5. 绘制给定系统的稳态响应

1. 给定序列，绘制其图像并分析其混叠

1.1 用MATLAB实现序列 x ( n ) = 3 c o s ( 0.125 π n + 0.2 π ) + 2 s i n ( 0.25 π n + 0.1 π ) x(n) = 3cos(0.125\pi n+0.2\pi)+2sin(0.25\pi n+0.1\pi) x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π)。

n=0:1:50; subplot(2,2,1);
x=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi);
stem(n,x,'fill'); xlabel('时间序列n '); title("原函数"); ylabel('x(n)'); axis([0 50 -5 4])

1.2 将1.1中的 x ( n ) x(n) x(n)扩展为以16为周期的函数 x ( n ) 16 = x ( n + 16 ) x(n)_{16}=x(n+16) x(n)16=x(n+16)，绘制出4个周期。

rem(a, b)函数的作用是求a/b的余数，本题要将序列扩展为16为周期的函数，只需要对自变量n进行取余操作，16为循环，便可以实现周期的扩展。

n=0:1:63; subplot(2,1,2);
x=3*cos(0.125*pi*rem(n,16)+0.2*pi)+2*sin(0.25*pi*rem(n,16)+0.1*pi);
stem(n,x,'fill'); xlabel ('时间序列n '); ylabel('x16(n)'); title("周期为16"); axis([0 64 -6 4])

1.3 将1.1中的 x ( n ) x(n) x(n)扩展为以10为周期的函数 x ( n ) 16 = x ( n + 10 ) x(n)_{16}=x(n+10) x(n)16=x(n+10)，绘制出4个周期。

n=0:1:39; subplot(2,2,2);
x=3*cos(0.125*pi*rem(n,10)+0.2*pi)+2*sin(0.25*pi*rem(n,10)+0.1*pi);
stem(n,x,'fill'); xlabel ('时间序列n '); ylabel('x10(n)'); title("周期为10");

2. 设计序列的移位和累加函数

x ( n ) = [ 1 , − 1 , 3 , 5 ] x(n)=[1,-1,3,5] x(n)=[1,−1,3,5]，产生并绘制出下列序列的样本：

2.1 x 1 ( n ) = 2 x ( n + 2 ) − x ( n − 1 ) − 2 x ( n ) x_1(n)=2x(n+2)-x(n-1)-2x(n) x1(n)=2x(n+2)−x(n−1)−2x(n)

2.2 x 2 ( n ) = n ∑ k = 1 5 x ( n − k ) x_2(n)=n\sum_{k=1}^5x(n-k) x2(n)=n∑k=15x(n−k)

序列移位就是 x ( n ) x(n) x(n)的每个样本都要移位一个量 k k k得到一个新的移位的序列 y ( n ) y(n) y(n)， y ( n ) = x ( n − k ) y(n)={x(n-k)} y(n)=x(n−k)

如果令 m = n − k m=n-k m=n−k，那么 n = m + k n = m + k n=m+k，上面运算给出为： y ( m + k ) = x ( m ) y(m+k)={x(m)} y(m+k)=x(m)，将此移位操作封装为函数：

function [y,n] = sigshift(x,m,k)
%implements y(n) = x(n - k)
%_________________________
%[y,n] = sigshift(x,m,k)
%
n = m+k;
y = x;

此外，我们需要将三个移位的子序列相加得到最终的序列，为此，建立函数 s i g a d d ( ) sigadd() sigadd()：

function [y,n] = sigadd(x1,n1,x2,n2)
% 实现序列相加
% 实现y(n)=x1(n)+x2(n)
% [y,n ]=sigadd(x1 ,n1,x2,n2)
n = min(min(n1),min(n2)):max(max(n1),max(n2)); % 求出序列的长度
y1 = zeros(1,length(n));
y2=y1;
y1 (find( (n>=min (n1) ) & (n<=max (n1) )==1) )=x1 ; % 将子序列的值复制到和序列的对应位置
y2 (find( (n>=min (n2) ) & (n<=max (n2) )==1) )=x2 ;
y=y1+y2;
end

最终调用以上两个子函数实现序列的相加并求和：

% x(n) = 2x(n+2)-x(n-1)-2x(n)
clc; clear;
n=0:3; x = [1 -1 3 5];
[x11, n11] = sigshift(x, n, 2);[x12, n12] = sigshift(x, n, -1);[x13, n13] = sigshift(x, n, 0);
[x1, n] = sigadd(2*x11, n11, (-1)*x12, n12);[x1, n] = sigadd(x1, n, (-2)*x13, n13);
subplot(2,1,1); stem(n, x1, 'fill');
xlabel ('时间序列n ');ylabel('x(n)');title("2x(n+2)-x(n-1)-2x(n)");
% x(n) = n(x(n-1)+x(n-2)+x(n-3)+x(n-4)+x(n-5))
n=0:3; x = n .* [1 -1 3 5];
[x11, n11] = sigshift(x, n, -1);[x12, n12] = sigshift(x, n, -2);[x13, n13] = sigshift(x, n, -3);
[x14, n14] = sigshift(x, n, -4);[x15, n15] = sigshift(x, n, -5);
[x1, n1] = sigadd(x11, n11, x12, n12);[x2, n2] = sigadd(x13, n13, x14, n14);
[x1, n1] = sigadd(x1, n1, x15, n15);[x1, n1] = sigadd(x1, n1, x2, n2);
subplot(2,1,2); stem(n1, x1 ,'fill');
xlabel ('时间序列n ');ylabel('x(n)');title("n(x(n-1)+x(n-2)+x(n-3)+x(n-4)+x(n-5))");

3. 绘制因果系统的幅频响应和相频响应

给定因果系统 H ( z ) = ( 1 + 2 z − 1 ) ( 1 − 0.67 z − 1 + 0.9 z − 2 ) H(z)=\frac{(1+\sqrt{2}z^{-1})}{(1-0.67z^{-1}+0.9z^{-2})} H(z)=(1−0.67z−1+0.9z−2)(1+2 z−1)，求出并绘制出其幅频响应和相频响应（提示使用MATLAB的freqz()计算频率响应）。

MATLAB画频率响应的函数为：

[h,w] = freqz(b,a,n)

b和a为传递函数系数，h为频率响应，w为角频率，范围 ( 0 , π ) (0, \pi) (0,π)，编程如下：

b = [2 0]; a = [1 -0.67 0.9];
[H,x] = freqz(b,a);
Hf=abs(H); % 取幅度值实部
Hx=angle(H); % 取相位值对应相位角
[AX, H1, H2] = plotyy(x, 20*log10(Hf), x, Hx); % 幅值变换为分贝单位
set(get(AX(1), 'Ylabel'), 'String', '幅频特性曲线');
set(get(AX(2), 'Ylabel'), 'String', '相频特性曲线');
title('因果系统的幅频特性和相频特性曲线');
set(H1, 'LineStyle', '--'); set(H2, 'LineStyle', ':');

4. 验证DTFT的折叠特性

令 x ( n ) x(n) x(n)为一复随机序列，其实部和虚部在 − 10 ≤ n ≤ 20 -10\leq n\leq 20 −10≤n≤20区间，均匀分布于[0,1]之间，用MATLAB验证DTFT的折叠特性（提示：序列的折叠用fliplr()函数，折叠后序列的下标为-fliplr(n)）。

由题，首先生成 − 10 ≤ n ≤ 20 -10\leq n\leq 20 −10≤n≤20共包含31个随机复数的序列：

% 折叠前的离散时间序列
n = -10:20; M = rand([1,31]); K = rand([1,31]); x = M + 1i*K;
k = 0:100; w = (pi/100)*k;

DTFT的折叠特性可以表示为： x ( − n ) = D T F T X ( e − j w ) x(-n)\stackrel{\mathrm{DTFT}}{=}X(e^{-jw}) x(−n)=DTFTX(e−jw)，大致的验证思路就是画出序列反褶后做傅里叶变换以及对原序列做傅里叶变换时将虚数部分取反，两种变换后的的幅度、相位、实数和虚数部分，比较是否相同。

%j->-j
X = x * (exp(1j * pi/100)).^(n' * k);
magX = abs(X); angX = angle(X); realX = real(X); imagX = imag(X);
subplot(4,2,1);plot(w/pi,magX,'r');title('幅度部分');xlabel('w/pi');ylabel('幅度');
subplot(4,2,3);plot(w/pi,angX,'r');title('相位部分');xlabel('w/pi');ylabel('相位');
subplot(4,2,5);plot(w/pi,realX,'r');title('实部');xlabel('w/pi');ylabel('实数值');
subplot(4,2,7);plot(w/pi,imagX,'r');title('虚部');xlabel('w/pi');ylabel('虚数值');
% 折叠后的离散时间序列
n = -fliplr(n); x = fliplr(x);
k = 0:100; w = (pi/100)*k;
X = x * (exp(-1j * pi/100)).^(n' * k);
magX = abs(X); angX = angle(X); realX = real(X); imagX = imag(X);
subplot(4,2,2);plot(w/pi,magX);title('幅度部分');xlabel('w/pi');ylabel('幅度');
subplot(4,2,4);plot(w/pi,angX);title('相位部分');xlabel('w/pi');ylabel('相位');
subplot(4,2,6);plot(w/pi,realX);title('实部');xlabel('w/pi');ylabel('实数值');
subplot(4,2,8);plot(w/pi,imagX);title('虚部');xlabel('w/pi');ylabel('虚数值');

5. 绘制给定系统的稳态响应

一个线性时不变系统的差分方程为 y ( n ) = 0.8 y ( n − 1 ) + x ( n ) y(n)=0.8y(n-1)+x(n) y(n)=0.8y(n−1)+x(n)，求出并画出它对输入 x ( n ) = c o s ( 0.05 π n ) u ( n ) x(n)=cos(0.05\pi n)u(n) x(n)=cos(0.05πn)u(n)的稳态响应。

由 y ( n ) − 0.8 y ( n − 1 ) = x ( n ) y(n)-0.8y(n-1)=x(n) y(n)−0.8y(n−1)=x(n)，得出频率响应 H ( e j w ) = 1 1 − 0.8 e − j w = e j w e j w − 0.8 H(e^{jw})=\frac{1}{1-0.8e^{-jw}}=\frac{e^{jw}}{e^{jw}-0.8} H(ejw)=1−0.8e−jw1=ejw−0.8ejw，

又 x ( n ) = c o s ( 0.05 π n ) , w = 0.05 , θ = 0 x(n)=cos(0.05\pi n), w=0.05,\theta =0 x(n)=cos(0.05πn),w=0.05,θ=0

H ( e j 0.05 π ) = 1 1 − 0.8 e − j 0.05 π ≈ 4.09 e − 0.54 j H(e^{j0.05\pi})=\frac{1}{1-0.8e^{-j0.05\pi}}\approx4.09e^{-0.54j} H(ej0.05π)=1−0.8e−j0.05π1≈4.09e−0.54j

% y(n) = 0.8y(n-1)+x(n), x(n) = cos(0.05πn)u(n)
clc; clear;
Y_vect = [1 -0.8]; X_vect = [1]; % 分子和分母的系数
% 绘制列向量Y的零点和列向量X的极点
subplot(3,2,1); zplane(Y_vect, X_vect);title("零点极点图");
% 时域脉冲响应图
subplot(3,2,2); impz(Y_vect, X_vect);title("脉冲响应图");
legend('h(n)');
t = 0 : 100;
Input = cos(0.05*pi*t); % 输入信号，频率w=0.05pi，相位为0
Output = filter(X_vect, Y_vect, Input); % 输入信号的响应
% x(n)
subplot(3,1,2);stem(Input);axis([0 100 -1 1]);title("x(n)原始输入");
% x(n)的响应图
subplot(3,1,3);stem(Output);axis([0 100 -5 5]); title("x(n)稳态响应图");