Python 3.11 零极点绘图实战：5种滤波器频率特性分析与代码复现

信号处理工程师在设计滤波器时，往往需要快速验证系统的频率响应特性。零极点分析作为频域设计的重要工具，能够直观反映滤波器的通带、阻带特性。本文将使用Python 3.11的科学计算栈（NumPy+SciPy+Matplotlib），通过5个典型滤波器案例，演示从零极点分布到频率响应可视化的完整工程实现。

1. 环境配置与基础理论

在开始前，确保已安装Python 3.11及以下库：

pip install numpy scipy matplotlib

零极点与系统函数的关系可通过传递函数表示：

H(s) = K * ∏(s-z_i)/∏(s-p_j)

其中 z_i 为零点， p_j 为极点， K 为增益系数。频率响应可通过 s=jω 代入计算。

关键工具函数 ：

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

def plot_response(z, p, k=1, fs=None):
    """绘制零极点图与频率响应"""
    plt.figure(figsize=(12, 6))
    
    # 零极点图
    plt.subplot(121)
    plt.scatter(np.real(z), np.imag(z), marker='o', facecolors='none', edgecolors='r', label='Zeros')
    plt.scatter(np.real(p), np.imag(p), marker='x', color='b', label='Poles')
    unit_circle = plt.Circle((0,0), 1, fill=False, linestyle='--', alpha=0.5)
    plt.gca().add_patch(unit_circle)
    plt.axhline(0, color='black', alpha=0.2)
    plt.axvline(0, color='black', alpha=0.2)
    plt.legend()
    
    # 频率响应
    plt.subplot(122)
    if fs:  # 离散系统
        w, h = signal.freqz_zpk(z, p, k, fs=fs)
        plt.plot(w, 20*np.log10(np.abs(h)), 'b')
    else:   # 连续系统
        w, h = signal.freqresp((z, p, k))
        plt.semilogx(w, 20*np.log10(np.abs(h)), 'b')
    plt.grid()
    return plt

2. 连续系统案例分析

2.1 二阶低通滤波器

特性 ：两个极点位于左半平面，无有限零点

poles = [-1 + 1j, -1 - 1j]  # 共轭极点
zeros = []                   # 无有限零点
k = 1                        # 增益

plt = plot_response(zeros, poles, k)
plt.title("2nd Order Lowpass")
plt.show()

输出特征 ：

• 幅频曲线在低频段平坦
• 高频段以-40dB/十倍频程衰减
• 相位从0°开始逐渐滞后

2.2 右半平面零点系统

特殊现象 ：零点位于右半平面导致非最小相位特性

zeros = [0.5 + 0j]  # 右半平面零点
poles = [-2 + 0j]   # 左半平面极点

plt = plot_response(zeros, poles)
plt.title("Non-minimum Phase System")
plt.show()

对比观察 ：

零点位置	相位特性	群延迟
左半平面	最小相位	较小
右半平面	非最小相位	较大

3. 离散系统案例分析

3.1 基本高通滤波器

Z域实现 ：零点在z=1，极点在原点

zeros = [1 + 0j]    # 单位圆上零点
poles = [0 + 0j]    # 原点极点
fs = 1000           # 采样率

plt = plot_response(zeros, poles, fs=fs)
plt.title("Discrete Highpass")
plt.show()

频响特点 ：

• 直流分量（ω=0）完全抑制
• 高频分量无衰减通过
• 相位响应呈线性变化

3.2 带通滤波器设计

设计要点 ：共轭极点靠近单位圆，零点合理配置

theta = np.pi/4  # 中心频率对应角度
r = 0.9          # 极点半径

zeros = [1, -1]  # 固定零点
poles = [r*np.exp(1j*theta), r*np.exp(-1j*theta)]  # 共轭极点

plt = plot_response(zeros, poles, fs=fs)
plt.title("Bandpass Filter")
plt.show()

参数优化建议 ：

1. 极点半径 r 越接近1，谐振峰越尖锐
2. 零点位置影响阻带衰减特性
3. 使用 signal.butter 可快速生成标准滤波器：

b, a = signal.butter(4, [0.2, 0.4], 'bandpass')

4. 综合实战：多级滤波器级联

设计一个满足以下指标的滤波器：

• 通带：0.1π ≤ ω ≤ 0.3π
• 阻带衰减 > 40dB
• 通带波纹 < 1dB

分步实现 ：

# 第一级：低通
z1, p1, k1 = signal.cheby1(4, 1, 0.3, 'lowpass', output='zpk')

# 第二级：高通 
z2, p2, k2 = signal.cheby1(4, 1, 0.1, 'highpass', output='zpk')

# 级联组合
z_total = np.concatenate([z1, z2])
p_total = np.concatenate([p1, p2])
k_total = k1 * k2

plt = plot_response(z_total, p_total, k_total, fs=fs)
plt.title("Cascaded Filter")
plt.show()

性能验证 ：

w, h = signal.freqz_zpk(z_total, p_total, k_total)
print(f"Passband ripple: {np.max(np.abs(h[10:30])) - np.min(np.abs(h[10:30])):.2f} dB")
print(f"Stopband attenuation: {-20*np.log10(np.max(np.abs(h[50:]))):.2f} dB")

5. 高级技巧与问题排查

5.1 数值稳定性处理

当极点非常接近单位圆时，可能出现计算溢出：

# 不稳定示例
poles = [0.9999*np.exp(1j*np.pi/3)]  
zeros = []

# 稳定化处理
poles = [0.99*np.exp(1j*np.pi/3)]  # 稍微远离单位圆

5.2 频率响应精细化控制

使用 worN 参数提高频率分辨率：

w, h = signal.freqz_zpk(z, p, k, worN=4096)

5.3 群延迟分析

w, gd = signal.group_delay((z, p, k))
plt.plot(w, gd)
plt.title('Group Delay')
plt.grid()

通过这5个案例的完整实现，我们建立了从理论分析到工程实践的完整链路。实际项目中，建议结合 scipy.signal 的现成滤波器设计函数，再通过零极点分析验证其特性