Python 信号与系统实战：6个零极点案例解析系统频率特性（附Bode图代码）

Python 信号与系统实战：6个零极点案例解析系统频率特性（附Bode图代码）

在信号处理领域，理解系统频率特性是每个工程师的必修课。传统教学中，我们往往通过纸笔计算和定性分析来掌握这一概念，但在实际工程应用中，如何快速验证理论分析、直观呈现系统行为才是真正的核心需求。本文将带你用Python这把"瑞士军刀"，从零极点分布出发，完整实现系统频率特性的自动化分析与可视化。

1. 理论基础与工具准备

零极点分析是理解线性时不变(LTI)系统频率响应的钥匙。系统函数H(s)或H(z)的分子分母多项式根（即零点和极点）在复平面的分布，直接决定了系统的幅频和相频特性。这种关系可以概括为：

• 极点 ：主导系统响应的峰值，越靠近虚轴（连续系统）或单位圆（离散系统），影响越显著
• 零点 ：形成频率响应的谷值，同样遵循距离原则
• 主零极点 ：距离关键区域（虚轴/单位圆）最近的零极点对系统行为起决定性作用

Python生态中， scipy.signal 模块提供了完整的信号处理工具链。我们需要重点关注以下函数：

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

# 核心函数
signal.ZerosPolesGain()  # 零极点增益模型
signal.bode()            # 连续系统波特图计算
signal.freqz()           # 离散系统频率响应

配置开发环境只需一行命令：

pip install numpy scipy matplotlib

2. 连续时间系统案例分析

2.1 典型带通滤波器（第一小题）

考虑系统函数： $$H(s) = \frac{s}{(s+2)(s+3)}$$

其零极点分布为：

• 零点：s = 0
• 极点：s = -2, s = -3

实现代码 ：

# 构建系统
num = [1, 0]       # s
den = [1, 5, 6]    # (s+2)(s+3) = s² + 5s + 6
sys = signal.TransferFunction(num, den)

# 计算频率响应
w, mag, phase = signal.bode(sys, np.logspace(-1, 2, 500))

# 绘制波特图
plt.figure(figsize=(10, 6))
plt.semilogx(w, mag, 'g', label='Magnitude (dB)')
plt.twinx()
plt.semilogx(w, phase, 'r', label='Phase (deg)')
plt.title('Bandpass Filter Frequency Response')
plt.grid(which='both', linestyle='--')
plt.legend()
plt.show()

关键观察 ：

1. 幅频特性呈现明显的钟形曲线，证实了带通特性
2. 相位从0°开始，随频率增加单调下降至-90°
3. -3dB带宽约为1.5 rad/s到4 rad/s

2.2 近极点主导系统（第二小题）

系统函数： $$H(s) = \frac{(s+100)}{(s+1)^2}$$

特性分析 ：

num = [1, 100]
den = [1, 2, 1]  # (s+1)² = s² + 2s + 1
sys = signal.TransferFunction(num, den)

w, mag, phase = signal.bode(sys, np.logspace(-2, 3, 500))

运行结果 ：

• 幅频曲线呈现典型低通特性
• 截止频率约在1 rad/s附近
• 远离虚轴的零点(s=-100)几乎不影响频率响应

提示：当零极点距离虚轴超过10倍以上时，其影响通常可以忽略不计

2.3 全通系统（第三小题）

特殊零极点配置： $$H(s) = \frac{(s-1)(s-2)}{(s+1)(s+2)}$$

代码实现要点 ：

zeros = [1, 2]    # 右半平面零点
poles = [-1, -2]  # 左半平面极点
k = 1             # 增益
sys = signal.ZerosPolesGain(zeros, poles, k)

特征验证 ：

1. 幅频响应平坦（全通特性）
2. 相位单调递减
3. 零极点关于虚轴对称配置是实现全通系统的充分条件

3. 离散时间系统案例分析

3.1 低通滤波器（第四小题）

离散系统函数： $$H(z) = \frac{z^2}{(z-0.9)}$$

数字信号处理要点 ：

• 单位圆上的零点不影响幅频响应
• 靠近单位圆的极点(0.9)产生低通特性

Python实现 ：

num = [1, 0, 0]  # z²
den = [1, -0.9]  
w, h = signal.freqz(num, den, worN=2000)

plt.plot(w, 20*np.log10(abs(h)))
plt.title('Discrete Lowpass Filter')
plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [rad/sample]')
plt.grid()

3.2 非最小相位系统（第五小题）

系统函数： $$H(z) = \frac{(z-1.2)^2}{(z-0.8e^{jπ/4})(z-0.8e^{-jπ/4})}$$

特殊性质 ：

• 单位圆外零点导致非最小相位
• 共轭极点产生谐振峰

关键代码 ：

zeros = [1.2, 1.2]  # 二阶零点
poles = [0.8*np.exp(1j*np.pi/4), 0.8*np.exp(-1j*np.pi/4)]
sys = signal.ZerosPolesGain(zeros, poles, 1)

w, h = signal.freqz(*sys.to_tf(), whole=True)

3.3 离散全通系统（第六小题）

零极点对称配置： $$H(z) = \frac{(z-1.5e^{jπ/6})(z-1.5e^{-jπ/6})}{(z-0.67e^{jπ/6})(z-0.67e^{-jπ/6})}$$

验证要点 ：

angle = np.pi/6
r = 0.67
zeros = [1.5*np.exp(1j*angle), 1.5*np.exp(-1j*angle)]
poles = [r*np.exp(1j*angle), r*np.exp(-1j*angle)]
sys = signal.ZerosPolesGain(zeros, poles, 1)

# 幅频响应验证
w, h = signal.freqz(*sys.to_tf())
assert np.allclose(abs(h), 1)  # 全通验证

4. 通用分析工具开发

基于上述案例，我们可以封装一个通用的频率特性分析函数：

def analyze_system(zeros, poles, gain, is_discrete=False, fs=1.0):
    """
    通用系统频率特性分析工具
    
    参数：
        zeros: 零点列表（复数）
        poles: 极点列表（复数）
        gain: 系统增益
        is_discrete: 是否为离散系统
        fs: 采样频率（离散系统需要）
    
    返回：
        fig: 可视化图形
        sys: 系统对象
    """
    sys = signal.ZerosPolesGain(zeros, poles, gain)
    fig = plt.figure(figsize=(12, 8))
    
    if is_discrete:
        # 离散系统分析
        w, h = signal.freqz(*sys.to_tf(), worN=2000, fs=fs)
        ax1 = fig.add_subplot(211)
        ax1.plot(w, 20*np.log10(abs(h)))
        ax1.set_ylabel('Magnitude [dB]')
        
        ax2 = fig.add_subplot(212)
        ax2.plot(w, np.angle(h))
        ax2.set_ylabel('Phase [rad]')
    else:
        # 连续系统分析
        w, mag, phase = signal.bode(sys)
        ax1 = fig.add_subplot(211)
        ax1.semilogx(w, mag)
        ax1.set_ylabel('Magnitude [dB]')
        
        ax2 = fig.add_subplot(212)
        ax2.semilogx(w, phase)
        ax2.set_ylabel('Phase [deg]')
    
    plt.xlabel('Frequency [rad/s]' if not is_discrete else 'Frequency [Hz]')
    plt.grid(True)
    return fig, sys

使用示例 ：

# 分析连续系统
zeros = [0]
poles = [-2, -3]
fig, _ = analyze_system(zeros, poles, 1)

# 分析离散系统
zeros = [1.2, 1.2]
poles = [0.8*np.exp(1j*np.pi/4), 0.8*np.exp(-1j*np.pi/4)]
fig, _ = analyze_system(zeros, poles, 1, is_discrete=True, fs=1000)

5. 工程实践中的注意事项

在实际应用中，有几个关键点需要特别注意：

1. 数值稳定性问题 ：

   • 高阶系统容易出现数值计算误差
   • 建议将系统转换为二阶节(SOS)形式：

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

2. 频率范围选择 ：

   • 连续系统：对数尺度通常更合适
   • 离散系统：关注0到π弧度/样本的范围

3. 可视化优化技巧 ：

   plt.figure(figsize=(10, 8))
   plt.subplot(2, 1, 1)
   plt.semilogx(w, mag, 'b', linewidth=2)
   plt.grid(which='both', alpha=0.5)
   plt.title('Bode Plot')
   

4. 常见问题排查表 ：

在完成多个项目后，发现最实用的调试方法是逐步验证：

1. 先确认零极点位置是否正确
2. 检查增益归一化
3. 验证频率范围设置是否合理
4. 最后处理可视化细节