信号与系统 MATLAB 代码

第 1 章 信号变换与运算

1.1 常用信号

​x
1
% 几种典型信号
2
% matlab 相比于 "傻瓜绘图软件",
3
% 优点: 代码可以复用, 函数种类众多.
4
% 缺点: 软件启动较慢, 操作相对繁琐.
5

6
% 1 指数
7
% 1.1 指数信号
8
t = -2 : 0.1 : 10;
9
ft = 2 * exp(-0.4*t);
10
plot(t, ft)
11
title('1.1 指数信号')
12

13
% 1.2 单边指数信号
14
t = -2 : 0.1 : 10;
15
ft = 2 * exp(-0.4*t) .* heaviside(t);
16
% 注意是点乘; heaviside(t) 是单位阶跃函数 (也可以当成示性函数使用)
17
% heaviside(t) 即 (t>0)
18
plot(t, ft)
19
title('1.2 单边指数信号')
20

21
% 1.3 复指数信号
22
t = 0 : 0.01 : 3;
23
z = exp((-1+10j)*t);
24
% 无定义时, i 和 j 都表示虚数单位, 并且与实数相乘无需写乘号
25
subplot 221, plot(t, real(z)), title('实部')
26
subplot 223, plot(t, imag(z)), title('虚部')
27
subplot 222, plot(t, abs(z)), title('模')
28
subplot 224, plot(t, angle(z)), title('相角')
29
% 第一次知道 subplot、plot 和 title 可以写在一行
30

31

32

33
% 2 正弦
34
% 2.1 正弦信号
35
t = 0 : 0.01 : 8;
36
ft = 1 * sin(2*pi*t + pi/6);
37
plot(t, ft)
38
title('2.1 正弦信号')
39

40
% 2.2 衰减正弦信号
41
t = 0 : 0.01 : 8;
42
ft = exp(-0.5*t) .* sin(2*pi*t + pi/6);
43
plot(t, ft)
44
title('2.2 衰减正弦信号')
45

46
% 2.3 抽样信号
47
t = -32 : 0.1 : 32;
48
ft = sin(t) ./ t;
49
plot(t, ft)
50
title('2.3 抽样信号')
51

52
% 2.4 辛格函数
53
t = -10 : 0.1 : 10;
54
ft = sin(pi*t) ./ (pi*t);
55
plot(t, ft)
56
title('2.4 辛格函数')
57

58

59

60
% 3 钟形信号
61
t = -5 : 0.1 : 5;
62
sigma = 1;
63
ft = exp(-t.^2 / (2*sigma^2)) / (sqrt(2*pi) * sigma);
64
plot(t, ft)
65
title('3 钟形信号')
66

67

68

69
% 4 奇异信号
70
% 4.1 单位斜边信号
71
t = -1 : 0.1 : 5;
72
ft = t .* heaviside(t);
73
plot(t, ft)
74
title('4.1 单位斜变信号')
75

76
% 4.2 单位阶跃信号
77
t = -1 : 0.001 : 3;  % 间隔小了, 精度会低
78
ft = heaviside(t);
79
plot(t, ft)
80
axis([-1, 3, -0.1, 1.2])
81
title('4.2 单位阶跃信号')
82

83
% 4.3 单位冲激信号
84
t1 = -1; t2 = 5;    % 信号起止时刻
85
t0 = 0;             % 信号平移量
86
dt = 0.01;          % 信号间隔时间
87
t = -1 : dt : 5;
88
x = zeros(1, length(t));
89
x(1, (-t0 - t1)/dt + 1) = 1 / dt;
90
stairs(t, x);
91
axis([t1, t2, 0, 1.2/dt])
92
title('4.3 单位冲激信号')
93

94
% 4.4 矩形脉冲信号
95
t = 0 : 0.001 : 4;
96
T = 1;
97
ft = rectpuls(t-2*T, 2*T);
98
% 其中 t-2*T 表示中心位置在 2*T
99
% 第二个参数表示矩形的宽度
100
plot(t, ft);
101
grid on
102
axis([0, 4, -0.5, 1.5])
103

104
% 4.5 周期性三角波脉冲信号
105
t = -5*pi : pi/10 : 5*pi;
106
ft = sawtooth(t, 0.3);
107
% 直译锯齿; 周期为 2*pi.
108
% 第二个参数 r 属于 [0, 1], 表示在 r * 2pi 处取得最大值.
109
plot(t, ft)
110
axis([-16, 16, -1.5, 1.5])
111
grid on

1.2 信号变换

xxxxxxxxxx
17
1
% 信号的运算
2
% 时移、反褶、尺度变换，可通过 y = subs(f, t, a*t - t_0) 实现
3
% 不建议使用 ezplot, 可替换为 fplot
4
% 如果不能使用 sym(), 可以直接输入引号内的符号表达式
5
syms t;
6
f1 = (4-t) * (heaviside(t) - heaviside(t-4));
7
f2 = sin(2*pi*t);
8
f3 = subs(f1, t, -t) + f1;
9
f4 = -f3;
10
f5 = f2 * f3;
11
f6 = f1 * f2;
12
subplot 321, fplot(f1), title('f1(t)')
13
subplot 322, fplot(f2), title('f2(t)')
14
subplot 323, fplot(f3), title('f1(-t) + f1(t)')
15
subplot 324, fplot(f4), title('-f1(-t) - f1(t)')
16
subplot 325, fplot(f5), title('f2(t) * f3(t)')
17
subplot 326, fplot(f6), title('f1(t) * f2(t)')

1.3 信号分解

xxxxxxxxxx
40
1
% 连续信号的分解
2
% 1 分解为直流分量和交流分量
3
% 符号函数不便于求均值
4
t = -2*pi : 0.01 : 2*pi;
5
ft = 1 + sin(t);
6
fd = mean(ft) * ones(1, length(t));
7
fa = ft - fd;
8
subplot 221, plot(t, ft), title('f(t)')
9
subplot 223, plot(t, fd), title('fd(t)')
10
subplot 224, plot(t, fa), title('fa(t)')
11
sgtitle('分解为直流分量和交流分量')
12

13
% 2 分解为奇分量和偶分量
14
syms t
15
ft = real((1 - (t - 1).^2) .^ 0.5);
16
fodd = real((ft - subs(ft, t, -t)) / 2);
17
feven = real((ft + subs(ft, t, -t)) / 2);
18
% 上述 real() 是必要的
19
subplot 221, fplot(ft), title('f(t)')
20
axis equal
21
axis([-2, 2, -1, 2])
22
subplot 223, fplot(fodd), title('fodd(t)')
23
axis equal
24
axis([-2, 2, -1, 2])
25
subplot 224, fplot(feven), title('feven(t)')
26
axis equal
27
axis([-2, 2, -1, 2])
28
% 受不了, 每次都要重新设置, 这么写代码太笨了
29
sgtitle('分解为奇分量和偶分量')
30

31
% 3 分解为脉冲分量
32
% 3.1 冲激信号
33
% 3.2 阶跃信号
34
% 这两个其实是把信号用积分的形式写了出来, 图像上没有差异
35

36
% 4 分解为实分量与虚分量
37
% 见 c1_signal.m 中的复指数信号部分, 那里已经画过了
38

39
% 5 分解为正交函数分量
40
% 由于会分解出无穷多个正交函数, 不方便全部展示, 也不画了

第 2 章 连续时域分析

2.1 系统响应

xxxxxxxxxx
12
1
% 连续时间系统的响应
2
% 以 2 y''(t) + y'(t) + 8 y(t) = f(t) 为例,
3
% 其中 f(t) 为冲激信号或阶跃信号.
4
% 只有单个元素的矩阵不建议使用方括号.
5
subplot 131, impulse(1, [2 1 8]), title('单位冲激信号')
6
subplot 132, step(1, [2 1 8]), title('单位阶跃信号')
7

8
% 对于 y''(t) + 2 y'(t) + y(t) = f'(t) + 2 f(t),
9
% 其中 f(t) = 5 e^(-2 t) u(t), 可以使用如下函数
10
t = 0 : 0.01 : 5;
11
f = 5 * exp(-2*t);
12
subplot 133, lsim([1 2], [1 2 1], f, t), title('一般信号')

2.2 函数卷积

同目录下 sconv.m 的内容:

xxxxxxxxxx
18
1
function [f, k] = sconv(f1, f2, k1, k2, p, isplot)
2
% 计算连续信号的卷积积分 f(t) = f1(t) * f2(t)
3
% k1 和 k2 分别是 f1 和 f2 的取样时刻
4
% p 为取样时间间隔
5
% plot 表示是否绘图
6
f = conv(f1, f2);
7
k0 = k1(1) + k2(1);
8
klen = length(k1) + length(k2) - 2;
9
k = k0 : p : (k0 + klen*p);
10

11
if nargin == 6 && isplot ~= 0
12
    subplot 221, plot(k1, f1), title('f1(t)')
13
    subplot 222, plot(k2, f2), title('f2(t)')
14
    subplot 223, plot(k, f), title('f(t) = f1(t) * f2(t)')
15
    h = get(gca, 'position');
16
    h(3) = 2.32 * h(3);
17
    set(gca, 'position', h)
18
end

具体的实例:

xxxxxxxxxx
27
1
% 连续时间信号的卷积积分
2
% 离散的可直接使用 conv(u, v)
3
% 连续的可取样后使用上述函数
4

5
% 例 2.16: t * I_{0 <= t <= 2} 与自身卷积
6
% 这里使用自己编写的 sconv 实现
7
p = 0.005;
8
k1 = 0 : p : 2;
9
f1 = 0.5 * k1;
10
[f, k] = sconv(f1, f1, k1, k1, p, true);
11

12
% 例 2.17
13
% 道理是一样的, 这里就不写了; 不过用 mathematica 写了.
14

15
% 例 2.18: I_{-1 <= t <= 1} 和 I_{0 <= t <= 1} 的卷积
16
% 这里不用 sconv, 而只使用其中部分代码
17
T = 0.01;   % 采样周期
18
n1 = -100 : 100;
19
f1 = ones(1, length(n1));
20
n2 = 0 : 100;
21
f2 = ones(1, length(n2));
22
y = T * conv(f1, f2);
23
n3 = n1(1) + n2(1) : n1(end) + n2(end);
24
plot(n3*T, y);
25
grid on;
26
xlabel('t/s');
27
ylabel('y(t)');

第 3 章 连续频域分析

3.1 周期信号的频谱

xxxxxxxxxx
34
1
% 3.1 周期信号的频谱
2

3
% 例 3.13: 周期矩形脉冲信号的双边频谱
4
syms t n        % T, A 和 tau 不需要用 syms 声明, 因为下面直接赋值了
5
T = 4; A = 1; tau = 1;
6

7
f = A * exp(-1j * 2*n*pi/T * t);    % 计算傅里叶级数时乘的函数
8
fn = int(f, t, -tau/2, tau/2) / T;  % 指数形式的傅里叶系数
9
% fn = simple(fn);                  % 书中说是用于化简
10
% 但是这个函数用不了, 甚至在官网都搜不到这个函数...不过也不需要用
11

12
n = [-20:-1, eps, 1:20];
13
evalfn = subs(fn, 'n', n);  % 书中参数写反了...
14

15
subplot 211, stem(n, evalfn, 'filled'), title('周期矩形脉冲的频谱');
16
subplot 212, stem(n, abs(evalfn), 'filled'), title('幅度谱');
17

18

19

20
% 例 3.14: 周期三角波信号的双边频谱
21
N = 10;
22
n1 = -N : -1;   % 小于 0 时的系数记为 c1
23
c1 = -4j * sin(n1*pi/2) / pi^2 ./ n1.^2;
24
c0 = 0;         % 等于 0 时的系数为 0
25
n2 = 1 : N;     % 大于 0 时的系数记为 c2
26
c2 = -4j * sin(n2*pi/2) / pi^2 ./ n2.^2;
27

28
cn = [c1, c0, c2];  % 傅里叶余弦级数系数
29
n = -N : N;
30
subplot 211, stem(n, abs(cn));
31
xlabel('\omega/\omega_0'); ylabel('C_n 的幅度');
32
subplot 212, stem(n, angle(cn));
33
xlabel('\omega/\omega_0'); ylabel('C_n 的相位');
34

3.2 非周期信号的频谱

xxxxxxxxxx
25
1
% 3.2 非周期信号的频谱
2

3
% 例 3.15: 单边指数信号的频谱 (符号计算)
4
syms t
5
f = exp(-2*t) * heaviside(t);
6
F = fourier(f);
7
subplot 211, fplot(abs(F)), grid on, title('幅频特性图')
8
subplot 212, fplot(angle(F)), grid on, title('相频特性图')
9

10
% 例 3.16: 门信号的频谱图 (数值计算)
11
R = 0.02;           % 采样间隔
12
t = -2 : R : 2;     % 采样时间点
13
N = 500;            % 样本容量
14
w = (0 : N)/N * (10*pi);    % 非负抽样点
15
W = [-fliplr(w), w(2:N+1)]; % 全部抽样点
16

17
f = heaviside(t+1) - heaviside(t-1);    % 门函数
18
F = f * exp(-1j*t'*w) * R;  % 非负抽样点
19
F = [fliplr(F), F(2:N+1)];  % 全部抽样点
20

21
subplot 211, plot(t, f, 'r');
22
xlabel('t'), ylabel('f(t)'), title('门函数 f(t)');
23
subplot 212, plot(W, real(F), 'b'); % 虽然 F 是实数, 但还是建议加上 real();
24
xlabel('w'), ylabel('F(w)'), title('傅里叶变换 F(w)');
25

3.3 傅里叶变换的性质

3.4 信号的采样过程

xxxxxxxxxx
21
1
% 3.4 信号的采样过程
2
% 例 3.19: 三角波的采样过程
3
syms t w f
4
f = (1 - 2*abs(t)) * exp(-1j * w * t);  % 被积函数
5
F = int(f, t, -1/2, 1/2);   % 指数形式的傅里叶系数
6

7
% 注意 fplot 第二个参数是二维向量, 而非采样点
8
subplot 311, fplot(abs(F), [-26*pi, 26*pi]), title('三角波的频谱');
9
axis([-26*pi, 26*pi, -0.1, 0.5]);
10
subplot 312, fplot(abs(F), [-4*pi, 4*pi]), title('低通滤波后的频谱');
11
axis([-26*pi, 26*pi, -0.1, 0.5]);
12

13
subplot 313, title('采样后的频谱'), hold on;
14
Ts = 0.2;       % 采样信号的周期
15
w0 = 2*pi / Ts; % 频谱的周期延拓
16
for k = -2 : 2
17
    ft = f * exp(-1j*w0*k*t);
18
    FT = int(ft, t, -1/2, 1/2);
19
    fplot(abs(FT)/Ts, [-4*pi-k*w0, 4*pi-k*w0]);
20
end
21

3.5 系统的频率特性

xxxxxxxxxx
24
1
% 3.5 系统的频率特性
2

3
% 例 3.20: 低通滤波器的频率响应
4
b = [0, 0, 1];
5
a = [0.08, 0.4, 1];
6
[h, w] = freqs(b, a, 100);
7
% 第三个参数为频率向量的频率点数 (默认为 200)
8
h1 = abs(h);
9
h2 = angle(h);
10
subplot 211, plot(w, h1), grid on;
11
xlabel('角频率 \omega'), ylabel('幅度'), title('幅频特性');
12
subplot 212, plot(w, h2*180/pi), grid on;
13
xlabel('角频率 \omega'), ylabel('相位 / 度'), title('相频特性');
14

15
% 例 3.21: 某连续 LTI 系统的微分方程
16
w = -2*pi : 0.1 : 2*pi;
17
b = [10, 5];
18
a = [1, 10, 8, 5];
19
H = freqs(b, a, w);
20
subplot 211, plot(w, abs(H)), title('幅频特性'), grid on;
21
xlabel('\omega(rad/s)'), ylabel('|H(j\omega)|');
22
subplot 212, plot(w, angle(H)), title('相频特性'), grid on;
23
xlabel('\omega(rad/s)'), ylabel('\phi(\omega)');
24

第 4 章 连续复频域分析

4.0 复频域常用函数

4.1 拉普拉斯变换

xxxxxxxxxx
39
1
% 4.1 拉普拉斯变换
2
% 例 4.27：两个信号的拉普拉斯变换
3
syms t w
4
% (1)
5
f1 = exp(-t) * cos(w*t);
6
F1 = laplace(f1);
7
% (2)
8
f2 = 3 * exp(-2*t);
9
F2 = laplace(f2);
10

11

12

13
% 例 4.28：绘制拉普拉斯变换的图形
14
syms t s
15
f = sin(t) * heaviside(t);
16
F = laplace(f);
17
% 得出 F = 1 / (s^2 + 1）
18

19
% 以下绘图代码可单独运行
20
% 绘制方法一：符号数学方法绘图
21
syms s x y
22
s = x + 1i*y;       % 虚数单位建议写为 1i 而非 i
23
Fs = 1 / (s^2+1);   % 注意这里要重新赋值为复变函数，因为要令 s = x + 1i * y
24
subplot 121, fmesh(abs(Fs));      % 网格曲面图；不建议使用 ezmesh()
25
subplot 122, fsurf(abs(Fs));      % 三维曲面图；不建议使用 ezsurf()
26
colormap(hsv);      % 设置图形中多条曲线的颜色顺序
27

28
% 绘制方法二：数值方法绘图
29
figure(2)
30
x1 = -5 : 0.1 : 5;
31
y1 = -5 : 0.1 : 5;
32
[x, y] = meshgrid(x1, y1);  % 产生网格矩阵 x, y
33
s = x + 1i*y;
34
Fs = 1 ./ (s.*s + 1);
35
subplot 121, mesh(x, y, abs(Fs));
36
subplot 122, surf(x, y, abs(Fs));
37
axis([-5, 5, -5, 5, 0, 8]);
38
colormap(hsv);
39

4.2 与傅氏变换关系

xxxxxxxxxx
50
1
% 4.2 傅氏变换与拉氏变换的关系
2
% 例 4.29 法一：符号绘图
3
clf;    % clear figure
4
syms x y s
5
s = x + 1i*y;
6
Fs = 1 ./ ((s+1)^2 + 1);
7

8
subplot 231, fmesh(abs(Fs));
9
view(-60, 20);  % 调整观察视角
10
axis([0, 5, -20, 20, 0, 0.5]);
11
title('拉氏变换网格图')
12

13
subplot 232, fsurf(abs(Fs));
14
view(-60, 20);  % 调整观察视角
15
axis([0, 5, -20, 20, 0, 0.5]);
16
title('拉氏变换曲面图')
17
colormap(hsv);
18

19
syms w
20
Fw = 1 ./ ((1i*w + 1)^2 + 1);
21
subplot 233, fplot(abs(Fw), [-20, 20]);
22
xlabel('频率 w')
23
title('傅氏变换幅频曲线')
24

25

26

27
% 例 4.29 法二：数值绘图
28
x1 = 0 : 0.1 : 5;
29
y1 = -20 : 0.1 : 20;
30
[x, y] = meshgrid(x1, y1);
31
s = x + 1i*y;
32
Fs = 1 ./ ((s+1).^2 + 1);
33

34
subplot 234, mesh(x, y, abs(Fs));
35
view(-60, 20);  % 调整观察视角
36
axis([0, 5, -20, 20, 0, 0.5]);
37
title('拉氏变换网格图')
38

39
subplot 235, surf(x, y, abs(Fs));
40
view(-60, 20);  % 调整观察视角
41
axis([0, 5, -20, 20, 0, 0.5]);
42
title('拉氏变换曲面图')
43
colormap(hsv);
44

45
w = -20 : 0.1 : 20;
46
Fw = 1 ./ ((1i*w + 1).^2 + 1);
47
subplot 236, plot(w, abs(Fw));
48
xlabel('频率 w')
49
title('傅氏变换幅频曲线')
50

4.3 拉普拉斯逆变换

xxxxxxxxxx
16
1
% 4.3 拉普拉斯逆变换
2
% 例 4.30：利用部分分式法进行拉普拉斯逆变换
3
a = conv([1, 0], [1, 1]);   % 分母的多项式系数
4
b = conv([1, 1], [1, 1]);   % 分子的多项式系数
5
den = conv(a, b);
6
[r, p] = residue([1, -2], den);
7

8
% 例 4.31：利用内置函数直接求拉普拉斯逆变换
9
syms s
10
% (1)
11
F1 = (2*s + 1) / (s^2 + 7*s + 10);
12
f1 = ilaplace(F1);
13
% (2)
14
F2 = s^2 / (s^2 + 3*s + 2);
15
f2 = ilaplace(F2);
16

4.4 系统函数零极点图

xxxxxxxxxx
35
1
% 4.4 系统函数的零、极点图
2
% 例 4.32：绘制零、极点图
3
b = [1, -1];
4
a = [1, 2, 2];
5
zs = roots(b);  % zeros
6
ps = roots(a);  % poles
7
plot(real(zs), imag(zs), 'o'); hold on;
8
plot(real(ps), imag(ps), 'rx', 'markersize', 12);
9
axis([-2, 2, -2, 2]); grid on;
10
legend('零点', '极点')
11

12
% 第五次翻转课堂练习题 2.6
13
z_real = [0, 1, 1];
14
z_imag = [0, -1, 1];
15
p_real = [-1, -1, 0, 0];
16
p_imag = [0, 0, -2, 2];
17
plot(z_real, z_imag, 'o'); hold on;
18
plot(p_real, p_imag, 'rx', 'markersize', 12);
19
axis([-3, 3, -3, 3]);
20
grid on;
21
legend('零点', '极点')
22

23
% 零极点的另一种画法：
24
H = tf([2 5 1], [1 3 5]);
25
pzmap(H)
26
grid on
27

28
pzmap(sys(:,:,1),'r',sys(:,:,2),'g',sys(:,:,3),'b')
29
sgrid
30

31
sys = tf([4.2,0.25,-0.004],[1,9.6,17]);
32
[p,z] = pzmap(sys);
33
plot(real(z), imag(z), 'o'); hold on;
34
plot(real(p), imag(p), 'rx', 'markersize', 12);
35

4.5 系统的复频域分析

xxxxxxxxxx
38
1
% 4.5 复频域分析
2
% 例 4.33：冲激响应的时域波形
3
% 备注：与微分算子、傅里叶变换那里的冲激响应使用函数相同
4
a = [1, 1, 16.25];
5
b = 1;  % 等价于 [1]
6
impulse(b, a, 5);
7
title('Impulse Response');
8
xlabel('Time');
9
ylabel('Amplitude');
10

11

12

13
% 例 4.34：零输入响应、零状态响应、全响应
14
% 首先对微分方程两端使用拉普拉斯变换，代入初值条件后化简，
15
syms t s
16
Yzi  = (3*s + 13) / (s^2 + 3*s + 2);    % 复频域 - 零输入响应
17
yzi = ilaplace(Yzi);                    % 时域 - 零输入响应
18

19
ft = 4 * exp(-2*t) * heaviside(t);      % 时域 - 激励信号
20
Fs = laplace(ft);                       % 复频域 - 激励信号
21

22
Yzs = Fs / (s^2 + 3*s + 2);             % 复频域 - 零状态响应
23
yzs = ilaplace(Yzs);                    % 时域 - 零状态响应
24

25
yt = simplify(yzi + yzs);               % 时域 - 全响应
26

27

28

29
% 例 4.35：给定电路，用拉普拉斯变换求电流响应
30
sys = tf(1, [0.1, 0.1]);    % 系统函数 H(s)
31
t = [0 : 0.01 : 10]';       % 抽样时间
32
e = sin(3*t);               % 激励信号
33
i = lsim(sys, e, t);        % 仿真电流信号
34
plot(t, e, '-.', t, i);     % 绘制激励和电流
35
xlabel('Time (sec)');
36
axis([0, 10, -4, 6]);
37
legend('e(t)', 'i(t)');
38

4.6 系统的频率响应图

xxxxxxxxxx
25
1
% 4.6 频率响应
2
% 例 4.36：系统函数的幅频特性和相频特性
3
b = 1;                  % 分子系数
4
a = [1, 0.4, 0.08];     % 分母系数
5
w = linspace(-pi, pi);  % 书中的 2*pi<512 疑似打错
6
h = freqs(b, a, w);         % 频率响应
7
subplot 211, plot(w, abs(h));    % 幅频特性曲线
8
subplot 212, plot(w, angle(h));  % 相频特性曲线
9

10
% 例 4.37：根据零极点分布画出幅频特性曲线
11
data = struct(...
12
    'title', {'(a)', '(b)', '(c)', '(d)'}, ...                  % 四个系统
13
    'zeros', {[], [0; 0], [-0.5], [1.2j; -1.2j]}, ...           % 零点坐标
14
    'poles', {[-2; -1], [-2, -1], [-2; -1], [-1+1j; -1-1j]}...  % 极点坐标
15
);  % 为绘图方便起见，使用结构体
16
omega = [0 : 0.01 : 6];     % 频率抽样点
17

18
for id = 1 : 4
19
    [b, a] = zp2tf(data(id).zeros, data(id).poles, 1);  % 零极点 -> 系统函数
20
    H = freqs(b, a, omega);                             % 指定频率点的响应
21
    subplot(2, 2, id), plot(omega, abs(H)), axis tight; % 绘制幅频特性曲线
22
    set(gca, 'YScale', 'log', 'FontSize', 16);          % 设置纵轴对数刻度
23
    title(data(id).title), xlabel('\omega'), ylabel('H(\omega)');
24
end
25

第 5 章 离散时域分析

5.1 离散时间信号的图像

xxxxxxxxxx
17
1
% 5.1 离散时间信号的图像
2
% 例 5.11：单边指数序列的波形图
3
n = 0 : 10;
4
a1 = 1.3;   x1 = a1.^n;
5
a2 = -1.3;  x2 = a2.^n;
6
a3 = 0.7;   x3 = a3.^n;
7
a4 = -0.7;  x4 = a4.^n;
8

9
subplot 221, stem(n, x1, 'fill'), grid on;
10
xlabel('n'), title('x(n) = 1.3^n');
11
subplot 222, stem(n, x2, 'fill'), grid on;
12
xlabel('n'), title('x(n) = -1.3^n');
13
subplot 223, stem(n, x3, 'fill'), grid on;
14
xlabel('n'), title('x(n) = 0.7^n');
15
subplot 224, stem(n, x4, 'fill'), grid on;
16
xlabel('n'), title('x(n) = -0.7^n');
17

5.2 离散时间信号的运算

xxxxxxxxxx
24
1
% 5.2 离散时间信号的运算
2
% 例 5.12：离散信号在基本运算后的波形图
3
% (1)   x1(n) = a^n [u(n) - u(n-N)]
4
a = 0.9;  N = 8;  n = -12 : 12;  n1 = n;
5
x1 = a.^n .* (n>=0 & n < N);    % n < N 不能取等，书中有误.
6
% 这样写是非常方便的, 也可以用 heavisdie(n) 表示 u(n)
7
subplot 411, stem(n1, x1, 'fill'), grid on;
8
title('x1(n)'), axis([-15 15 0 1]);
9

10
% (2)   x2(n) = x1(n + 3);
11
n2 = n1 - 3;  x2 = x1;  % 注意这里直接写为 n1 - 3
12
subplot 412, stem(n2, x2, 'fill'), grid on;
13
title('x2(n)'), axis([-15 15 0 1]);
14

15
% (3)   x3(n) = x1(n - 2);
16
n3 = n1 + 2;  x3 = x1;
17
subplot 413, stem(n3, x3, 'fill'), grid on;
18
title('x3(n)'), axis([-15 15 0 1]);
19

20
% (4)   x4(n) = x1(-n);
21
n4 = n1;   x4 = fliplr(x1);     % 这里书中有误, n4 与 x4 翻一次就可以了.
22
subplot 414, stem(n4, x4, 'fill'), grid on;
23
title('x4(n)'), axis([-15 15 0 1]);
24

5.3 系统的单位样值响应

xxxxxxxxxx
7
1
% 5.3 系统的单位样值响应
2
% 例 5.13：差分方程的单位样值响应
3
a = [3 -4 2];   % Y(n) 的系数
4
b = [1 2];      % X(n) 的系数
5
n = 0 : 50;     % 计算 51 个点
6
impz(b, a, 50); % 绘出单位样值响应
7

5.4 离散时间序列的卷积

xxxxxxxxxx
28
1
% 5.4 离散时间序列的卷积
2
% 例 5.14：离散序列的卷积和及其波形图
3
% x(n) = u(n) - u(n-10);
4
% h(n) = u(n) - u(n-5);
5
n = -10 : 10;   % 无需使用中括号
6
x = zeros(1, length(n));
7
x(n>=0 & n<10) = 1;     % 可以直接使用逻辑序列，而无需使用 find 函数
8
                        % 即时使用 find 函数，也无需在两侧加中括号
9
h = zeros(1, length(n));
10
h(n>=0 & n<5) = 1;
11

12
n1 = fliplr(-n);    h1 = fliplr(h);                     % 过程 1
13
n2 = n1;    h2 = [0, h1];   h2(length(h2)) = [];        % 过程 2
14
n3 = n2;    h3 = [0, h2];   h3(length(h3)) = [];        % 过程 3
15
n4 = -n;    nmin = min(n1) - max(n4);
16
nmax = max(n1) - min(n4);   ny = nmin : nmax;
17
y = conv(x, h);     % 卷积结果
18

19
subplot 321, stem(n, x, '*k'), title('x(n)');
20
subplot 322, stem(n, h, 'ok'), title('h(n)');
21
subplot 323, stem(n, x, '*k'), hold on;
22
             stem(n1, h1, 'k'), title('p1');
23
subplot 324, stem(n, x, '*k'), hold on;
24
             stem(n2, h2, 'k'), title('p2');
25
subplot 325, stem(n, x, '*k'), hold on;
26
             stem(n3, h3, 'k'), title('p3');
27
subplot 326, stem(ny, y, '.k'), title('y(n)');
28

5.5 离散时间系统的响应

xxxxxxxxxx
24
1
% 5.5 离散时间系统的响应
2
% 例 5.15：输入序列的时域波形图与 LTI 系统的零状态响应
3
a = [3 -4 2];   % y(n-k) 的系数
4
b = [1 2];      % x(n-k) 的系数
5
n = 0 : 20;
6
x = 0.5 .^ n;
7
y = filter(b, a, x);
8
subplot 211, stem(n, x, 'fill'), title('Input Seqence');
9
xlabel('n'); ylabel('x(n)');
10
subplot 212, stem(n, y, 'fill'), title('Response Sequence');
11
xlabel('n'); ylabel('y(n)');
12

13
% 例 5.16：系统的零状态响应
14
k1 = 0 : 7;     f1 = 0.8 .^ k1 .* ones(1, 8);
15
k2 = 0 : 3;     f2 = ones(1, 4);
16
% [f, k] = dconv(f1, f2, k1, k2);
17
f = conv(f1, f2);   % 教材中给的代码跑不了，换成这个了
18

19
% 教材中没有给出绘图代码. 这里的绘图功能是第一次遇到的
20
% 其中第三个图像, 可以省略横轴的 k3.
21
subplot 221, stem(k1, f1, '*b'), title('f1(n)');
22
subplot 222, stem(k2, f2, 'ob'), title('f2(n)');
23
subplot(2, 2, [3 4]), stem(f, 'ob'), title('f(n)');
24

第 6 章 离散 z 域分析

6.1 离散系统零极点图

6.2 系统单位样值响应

6.3 离散系统差分方程

6.4 z 域部分分式展开

6.5 z 变换与 z 逆变换