Python 3.11 信号分析实战：6个系统零极点分布与频率特性可视化（附代码）

在数字信号处理领域，理解系统的零极点分布与频率响应特性是掌握滤波器设计和系统分析的关键。本文将带你通过Python 3.11环境下的SciPy和Matplotlib工具包，将抽象的零极点理论转化为直观的可视化结果，特别适合正在学习《信号与系统》课程或需要快速验证系统特性的工程师。

1. 理论基础与工具准备

零极点分析是研究线性时不变(LTI)系统频率响应的核心方法。系统函数的零点和极点在复平面上的位置，直接决定了系统的幅频特性和相频特性。极点会使频率响应在该频率附近出现峰值，而零点则会导致频率响应在该频率附近出现谷值。

必备工具安装 ：

pip install numpy scipy matplotlib control

核心工具包功能对比 ：

工具包	主要功能	关键API
NumPy	基础数值计算	np.poly, np.roots
SciPy	科学计算	scipy.signal.bode
Matplotlib	数据可视化	plt.subplot, plt.semilogx
Control	控制系统分析	control.tf, control.pzmap

提示：本文所有代码示例均基于Python 3.11验证通过，建议使用Jupyter Notebook进行交互式实验

2. 连续时间系统分析

2.1 典型带通滤波器特性

我们先分析一个具有两个极点和一个零点的二阶系统，其系统函数为： $$H(s) = \frac{s}{(s+2)(s+3)}$$

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 定义系统
num = [1, 0]  # s
den = np.poly([-2, -3])  # (s+2)(s+3)
sys = signal.TransferFunction(num, den)

# 绘制零极点图
plt.figure(figsize=(10, 4))
plt.subplot(121)
zeros = np.roots(num)
poles = np.roots(den)
plt.scatter(np.real(zeros), np.imag(zeros), marker='o', facecolors='none', edgecolors='r', s=100)
plt.scatter(np.real(poles), np.imag(poles), marker='x', color='r', s=100)
plt.axhline(0, color='k', linestyle='--', alpha=0.3)
plt.axvline(0, color='k', linestyle='--', alpha=0.3)
plt.title('Pole-Zero Plot')
plt.xlabel('Real')
plt.ylabel('Imaginary')
plt.grid(True)

# 绘制波特图
plt.subplot(122)
w, mag, phase = signal.bode(sys)
plt.semilogx(w, mag, label='Magnitude')
plt.semilogx(w, phase, label='Phase')
plt.title('Bode Plot')
plt.xlabel('Frequency [rad/s]')
plt.ylabel('Magnitude [dB] / Phase [deg]')
plt.legend()
plt.grid(which='both', linestyle='--')
plt.tight_layout()
plt.show()

关键观察点 ：

• 零点位于原点，极点在负实轴(-2, -3)
• 幅频特性呈现典型的带通特征
• 相位从+90°开始单调下降至-90°

2.2 主极点主导的低通系统

当系统存在远离虚轴的零极点时，系统的频率特性主要由靠近虚轴的零极点决定。例如系统函数： $$H(s) = \frac{(s+100)}{(s+1)^2}$$

num = [1, 100]
den = [1, 2, 1]
sys = signal.TransferFunction(num, den)

# 可视化代码与上例类似...

特性分析 ：

• 极点-1比零点-100更靠近虚轴
• 幅频特性呈现二阶低通特征
• 高频段幅值下降速率为-40dB/decade

3. 离散时间系统分析

3.1 低通数字滤波器

离散系统分析使用z变换，零极点位于单位圆内表示稳定系统。分析系统： $$H(z) = \frac{z^2}{(z-0.9)}$$

# 离散系统定义
numz = [1, 0, 0]  # z^2
denz = [1, -0.9]  # (z-0.9)
dt = 0.1  # 采样时间
sysd = signal.TransferFunction(numz, denz, dt=dt)

# 频率响应计算
w, h = signal.freqz(numz, denz)
freq = w / (2*np.pi*dt)

plt.figure(figsize=(10, 4))
plt.subplot(121)
# 零极点图绘制...
plt.subplot(122)
plt.plot(freq, 20*np.log10(np.abs(h)), label='Magnitude')
plt.plot(freq, np.angle(h), label='Phase')
plt.title('Frequency Response')
plt.xlabel('Normalized Frequency')
plt.ylabel('Magnitude [dB] / Phase [rad]')
plt.legend()
plt.grid()
plt.show()

特性说明 ：

• 极点在0.9处靠近单位圆
• 零点在原点(二阶)
• 幅频特性呈现低通特征
• 相位呈现非线性变化

3.2 全通系统特性

当零极点关于单位圆呈共轭对称分布时，系统呈现全通特性。例如： $$H(z) = \frac{z^{-1}-a^*}{1-az^{-1}}, \quad |a|<1$$

a = 0.8*np.exp(1j*np.pi/4)  # 复数极点
numz = [1, -np.conj(a)]      # 零点在1/a*
denz = [1, -a]               # 极点在a

全通系统特征 ：

• 幅频响应为常数
• 相位响应非线性变化
• 常用于相位均衡

4. 综合案例：滤波器设计验证

结合零极点配置设计一个带阻滤波器，中心频率1kHz，采样率10kHz：

fs = 10000
f0 = 1000
w0 = 2*np.pi*f0/fs

# 在单位圆上放置零点，附近放置极点
zeros = [np.exp(1j*w0), np.exp(-1j*w0)]
poles = 0.95 * np.array(zeros)  # 极点在半径0.95处

num = np.poly(zeros)
den = np.poly(poles)

# 频率响应可视化
w, h = signal.freqz(num, den, worN=2000)
freq = w * fs / (2*np.pi)
plt.plot(freq, 20*np.log10(np.abs(h)))
plt.axvline(f0, color='r', linestyle='--')
plt.title('Bandstop Filter Response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude [dB]')
plt.grid()

5. 高级可视化技巧

5.1 三维频率响应展示

from mpl_toolkits.mplot3d import Axes3D

# 创建网格
w = np.linspace(0, np.pi, 500)
z = np.exp(1j*w)
H = np.polyval(num, z) / np.polyval(den, z)

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(121, projection='3d')
ax.plot_surface(np.real(z), np.imag(z), np.abs(H), cmap='viridis')
ax.set_title('Magnitude Response')
ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')

ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(np.real(z), np.imag(z), np.angle(H), cmap='coolwarm')
ax.set_title('Phase Response')
plt.tight_layout()

5.2 交互式零极点编辑器

from matplotlib.widgets import Slider

fig, ax = plt.subplots(figsize=(8, 6))
plt.subplots_adjust(bottom=0.25)

# 初始零极点
initial_pole = 0.8*np.exp(1j*np.pi/4)
pole_marker, = ax.plot(np.real(initial_pole), np.imag(initial_pole), 'rx', markersize=10)
ax.add_patch(plt.Circle((0,0), 1, color='b', fill=False))
ax.axis('equal')
ax.grid()

# 创建滑动条
axcolor = 'lightgoldenrodyellow'
ax_real = plt.axes([0.25, 0.15, 0.65, 0.03], facecolor=axcolor)
ax_imag = plt.axes([0.25, 0.1, 0.65, 0.03], facecolor=axcolor)

s_real = Slider(ax_real, 'Real', -1, 1, valinit=np.real(initial_pole))
s_imag = Slider(ax_imag, 'Imag', -1, 1, valinit=np.imag(initial_pole))

def update(val):
    pole = s_real.val + 1j*s_imag.val
    pole_marker.set_data([np.real(pole)], [np.imag(pole)])
    fig.canvas.draw_idle()

s_real.on_changed(update)
s_imag.on_changed(update)

6. 实际工程注意事项

1. 数值稳定性问题 ：

   • 高阶系统建议使用二阶节(SOS)形式

   sos = signal.zpk2sos(zeros, poles, gain)
   

2. 频率响应计算优化 ：

   # 使用对数间隔频率点
   w = np.logspace(-2, 2, 500)
   

3. 系统实现验证 ：

   # 时域响应验证
   t, y = signal.step(sys)
   plt.plot(t, y)
   

4. 滤波器设计实用函数 ：

   # 直接设计Butterworth滤波器
   order = 4
   cutoff = 0.2
   b, a = signal.butter(order, cutoff)
   

通过这6个完整的案例，我们展示了如何从零极点分布预测系统特性，并通过Python代码实现可视化验证。这种"理论分析→代码实现→可视化验证"的工作流程，正是现代信号处理工程师的必备技能。