由系统函数求零极点、频率响应(幅频特性、相频特性)的 Matlab 和 Python 方法
Author: Sijin Yu
本文以离散信号为例.
1. Matlab
1.1 tf2zpk() 函数
使用 tf2zpk() 函数可以获得频率响应的零极点.
matlab 的官方文档对 tf2zpk() 函数的用法介绍如下.
即, 给定系统函数
H
(
z
)
=
b
0
+
b
1
z
−
1
+
⋯
+
b
n
z
−
n
a
0
+
a
1
z
−
1
+
⋯
+
a
n
z
−
m
=
∑
i
=
0
n
b
i
z
−
i
∑
i
=
0
m
a
i
z
−
i
,
H(z)=\frac{b_0+b_1z^{-1}+\cdots+b_nz^{-n}}{a_0+a_1z^{-1}+\cdots+a_nz^{-m}}=\frac{\sum^{n}_{i=0}b_iz^{-i}}{\sum^m_{i=0}a_iz^{-i}},
H(z)=a0+a1z−1+⋯+anz−mb0+b1z−1+⋯+bnz−n=∑i=0maiz−i∑i=0nbiz−i,
令序列
b
=
[
b
0
,
b
1
,
⋯
,
b
n
]
,
a
=
[
a
0
,
a
1
,
⋯
,
a
m
]
b=[b_0,b_1,\cdots,b_n],a=[a_0,a_1,\cdots,a_m]
b=[b0,b1,⋯,bn],a=[a0,a1,⋯,am]. 系统函数的零极点表达式为
H
(
z
)
=
k
(
z
−
z
1
)
(
z
−
z
2
)
⋯
(
z
−
z
N
)
(
z
−
p
1
)
(
z
−
p
2
)
⋯
(
z
−
p
M
)
=
k
∏
i
=
1
N
(
z
−
z
i
)
∏
i
=
1
M
(
z
−
p
i
)
,
H(z)=k\frac{(z-z_1)(z-z_2)\cdots(z-z_N)}{(z-p_1)(z-p_2)\cdots(z-p_M)}=k\frac{\prod^N_{i=1}(z-z_i)}{\prod^M_{i=1}(z-p_i)},
H(z)=k(z−p1)(z−p2)⋯(z−pM)(z−z1)(z−z2)⋯(z−zN)=k∏i=1M(z−pi)∏i=1N(z−zi),
序列
z
=
[
z
1
,
z
2
,
⋯
,
z
N
]
,
p
=
[
p
1
,
p
2
,
⋯
,
p
M
]
z=[z_1,z_2,\cdots,z_N],p=[p_1,p_2,\cdots,p_M]
z=[z1,z2,⋯,zN],p=[p1,p2,⋯,pM] 分别表示
H
(
z
)
H(z)
H(z) 的零点和极点. 函数 [z, p, k]=tf2zpk(b, a) 返回以 b 和 a 序列为参数的系统方程的零点序列 z、极点序列 p、增益 k.
1.2 zplane() 函数
matlab 的官方文档对 zplane() 函数的用法介绍如下.
zplane() 有两个主要用法:
-
zplane(z, p). 传入参数为零点序列z和极点序列p. 直接作零极点图. -
zplane(b, a). 传入参数为系统函数的参数b和a. 函数自动根据系统函数作零极点图.\
在下文我们会验证这两种做法是等效的.
1.3 freqz() 函数
matlab 的官方文档对 freqz() 函数的用法介绍如下.
即, 给定系统函数
H
(
z
)
=
b
0
+
b
1
z
−
1
+
⋯
+
b
n
z
−
n
a
0
+
a
1
z
−
1
+
⋯
+
a
n
z
−
m
=
∑
i
=
0
n
b
i
z
−
i
∑
i
=
0
m
a
i
z
−
i
,
H(z)=\frac{b_0+b_1z^{-1}+\cdots+b_nz^{-n}}{a_0+a_1z^{-1}+\cdots+a_nz^{-m}}=\frac{\sum^{n}_{i=0}b_iz^{-i}}{\sum^m_{i=0}a_iz^{-i}},
H(z)=a0+a1z−1+⋯+anz−mb0+b1z−1+⋯+bnz−n=∑i=0maiz−i∑i=0nbiz−i,
令序列
b
=
[
b
0
,
b
1
,
⋯
,
b
n
]
,
a
=
[
a
0
,
a
1
,
⋯
,
a
m
]
b=[b_0,b_1,\cdots,b_n],a=[a_0,a_1,\cdots,a_m]
b=[b0,b1,⋯,bn],a=[a0,a1,⋯,am]. 函数 [h, w]=freqz(b, a) 返回频率响应序列 h 和频率序列 w. h 为复序列, 模 abs(h) 为幅频响应序列, 角 angle(h) 为相频响应序列.
1.4 Example
作下面系统函数的零极点图、幅频特性曲线、相频特性曲线.
H
(
z
)
=
0.5
z
−
1
−
0.72
z
−
3
3
+
2
z
−
1
−
0.87
z
−
2
.
H(z)=\frac{0.5z^{-1}-0.72z^{-3}}{3+2z^{-1}-0.87z^{-2}}.
H(z)=3+2z−1−0.87z−20.5z−1−0.72z−3.
代码如下:
% test.m
% author: Sijin Yu
clear;
figure(1);
b = [0 0.5 0 -0.72];
a = [3 2 0.87];
[z, p, k] = tf2zpk(b, a);
% -----零极点-----
subplot(2, 2, 1);
zplane(z, p);
title('zplane(z, p)');
subplot(2, 2, 2);
zplane(b, a);
title('zplane(b, a)');
% -----幅频特性-----
[h, w] = freqz(b, a);
H_abs = abs(h);
subplot(2, 2, 3);
plot(w, 20 * log10(H_abs)); % 以分贝为单位
title('幅频特性');
% -----相频特性-----
H_angle = angle(h);
subplot(2, 2, 4);
plot(w, H_angle);
title('相频特性');
结果如下:
2. Python
2.1 scipy.signal.tf2zpk() 函数
该函数依赖 scipy 库, 使用前应执行 pip install scipy.scipy 的官方文档介绍如下.
其用法和 matlab 中的 tf2zpk() 用法非常相似, 不多赘述.
2.2 zplane() 函数的自定义
Python 中没有直接实现 matlab 中 zplane() 函数的功能. (至少我没找到, 有知道的大佬欢迎留言.) 因此我自己实现了 zplane() 函数.
函数定义的代码如下
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Circle
def zplane(z, p, fig=None, ax=None):
if fig==None or ax==None:
fig, ax = plt.subplots(figsize=(4, 4))
circle = Circle(xy = (0.0, 0.0), radius = 1, alpha = 0.9, facecolor = 'white')
# 作单位园
theta = np.linspace(0, 2 * np.pi, 200)
x = np.cos(theta)
y = np.sin(theta)
ax.add_patch(circle)
ax.plot(x, y, color="darkred", linewidth=2)
lim = max(max(z), max(p), 1) + 1
# 控制坐标轴范围
plt.xlim([-lim, lim])
plt.ylim([-lim, lim])
# 作零极点
for i in z:
ax.plot(np.real(i),np.imag(i), 'bo')
for i in p:
ax.plot(np.real(i),np.imag(i), 'bx')
2.3 scipy.signal.freqz() 函数
scipy 库官方文档对其描述如下.
(由于内容过长, 不截图)
scipy.signal.freqz(b, a=1, worN=512, whole=False, plot=None, fs=6.283185307179586, include_nyquist=False)
Parameters:
-
b:array_like
Numerator of a linear filter. Ifbhas dimension greater than 1, it is assumed that the coefficients are stored in the first dimension, andb.shape[1:],a.shape[1:], and the shape of the frequencies array must be compatible for broadcasting. -
a:array_like
Denominator of a linear filter. Ifbhas dimension greater than 1, it is assumed that the coefficients are stored in the first dimension, andb.shape[1:],a.shape[1:], and the shape of the frequencies array must be compatible for broadcasting. -
worN:{None, int, array_like}, optional
If a single integer, then compute at that many frequencies (default is N=512). This is a convenient alternative to:np.linspace(0, fs if whole else fs/2, N, endpoint=include_nyquist).
Using a number that is fast for FFT computations can result in faster computations (see Notes).
If anarray_like, compute the response at the frequencies given. These are in the same units asfs. -
whole:bool, optional
Normally, frequencies are computed from 0 to the Nyquist frequency,fs/2(upper-half of unit-circle). Ifwholeis True, compute frequencies from 0 tofs. Ignored if worN isarray_like. -
plot:callable
Acallablethat takes two arguments. If given, the return parameterswandhare passed to plot. Useful for plotting the frequency response insidefreqz. -
fs:float, optional
The sampling frequency of the digital system. Defaults to 2*pi radians/sample (so w is from 0 to pi). -
include_nyquist:bool, optional
Ifwholeis False andworNis an integer, setting include_nyquist to True will include the last frequency (Nyquist frequency) and is otherwise ignored.
Returns:
-
w:ndarray
The frequencies at whichhwas computed, in the same units asfs. By default,wis normalized to the range [0, pi) (radians/sample). -
h:ndarray
The frequency response, as complex numbers.
2.4 Example
我们用回 1.4 中的例子.
代码如下:
# test.py
# author: Sijin Yu
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from matplotlib.patches import Circle
def zplane(z, p, fig=None, ax=None):
if fig==None or ax==None:
fig, ax = plt.subplots(figsize=(4, 4))
circle = Circle(xy = (0.0, 0.0), radius = 1, alpha = 0.9, facecolor = 'white')
# 作单位园
theta = np.linspace(0, 2 * np.pi, 200)
x = np.cos(theta)
y = np.sin(theta)
ax.add_patch(circle)
ax.plot(x, y, color="darkred", linewidth=2)
lim = max(max(z), max(p), 1) + 1
# 控制坐标轴范围
plt.xlim([-lim, lim])
plt.ylim([-lim, lim])
# 作零极点
for i in z:
ax.plot(np.real(i),np.imag(i), 'bo')
for i in p:
ax.plot(np.real(i),np.imag(i), 'bx')
b = np.array([0, 0.5, 0, -0.72])
a = np.array([3, 2, 0.87])
# -----零极点图-----
z, p, k = signal.tf2zpk(b, a)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(2, 2, 1)
zplane(z, p, fig=fig, ax=ax)
# -----幅频特性-----
w, h = signal.freqz(b, a)
ax = fig.add_subplot(2, 2, 3)
ax.plot(w, 20 * np.log10(abs(h))) # 以分贝为单位
# -----相频特性-----
ax = fig.add_subplot(2, 2, 4)
ax.plot(w, np.angle(h))
fig.savefig("result_py.jpg")
结果如下:
文章来源：https://www.toymoban.com/news/detail-489948.html
3. 总结
对比 matlab 的图和 Python 的图, 发现为原点的极点 Python 没有画出, 而 matlab 画出. 其余结果 matlab 与 Python 一致.文章来源地址https://www.toymoban.com/news/detail-489948.html
到了这里，关于由系统函数求零极点图、频率响应(幅频特性、相频特性)的 Matlab 和 Python 方法的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
