信号处理抽取多项滤波的数学推导与仿真

昨天的《 》，已经介绍了插值抽取的多项滤率，今天详细介绍多项滤波的数学推导，并附上实战仿真代码。

一、数学变换推导

1 多相分解的核心思想

将FIR滤波器的系数 h ( n ) h(n) h(n)按相位分组，每组对应输入信号的不同抽样相位。通过分相、滤波、重组，实现与原FIR等效的处理。

2 数学变换推导

FIR滤波器的系统函数可表示为： H ( z ) = ∑ n = 0 N − 1 h ( n ) z − n H(z) = \sum_{n=0}^{N-1} h(n) z^{-n} H(z)=n=0∑N−1​h(n)z−n 其中， h ( n ) h(n) h(n)为滤波器系数， N N N为阶数。 设分解因子为 M M M，则第 k k k个子滤波器系数为： h k ( m ) = h ( k + m M ) , 0 ≤ k < M h_k(m) = h(k + mM), \quad 0 \leq k < M hk​(m)=h(k+mM),0≤k<M 将FIR滤波器拆分为 M M M个并行的子滤波器（即多相分量 ），每个子滤波器处理信号的特定相位分量，再通过延迟和组合实现等效。目标形式为： H ( z ) = ∑ k = 0 M − 1 z − k H k ( z M ) H(z) = \sum_{k=0}^{M-1} z^{-k} H_k(z^M) H(z)=k=0∑M−1​z−kHk​(zM) 其中， H k ( z M ) H_k(z^M) Hk​(zM)表示第 k k k个子滤波器的系统函数。 将 H ( z ) H(z) H(z)按 M M M的整数倍延迟展开： H ( z ) = ∑ n = 0 N − 1 h ( n ) z − n = ∑ k = 0 M − 1 ∑ m = 0 K − 1 h ( k + m M ) z − ( k + m M ) H(z) = \sum_{n=0}^{N-1} h(n) z^{-n} = \sum_{k=0}^{M-1} \sum_{m=0}^{K-1} h(k + mM) z^{-(k + mM)} H(z)=n=0∑N−1​h(n)z−n=k=0∑M−1​m=0∑K−1​h(k+mM)z−(k+mM) 其中， K = ⌈ N M ⌉ K = \lceil \frac{N}{M} \rceil K=⌈MN​⌉（向上取整）。 将原系数按 M M M个相位分组：

• 第 k k k个子滤波器的系数为： h k ( m ) = h ( k + m M ) h_k(m) = h(k + mM) hk​(m)=h(k+mM)
• 其系统函数为： H k ( z M ) = ∑ m = 0 K − 1 h k ( m ) z − m M H_k(z^M) = \sum_{m=0}^{K-1} h_k(m) z^{-mM} Hk​(zM)=m=0∑K−1​hk​(m)z−mM 将 H ( z ) H(z) H(z)重写为： H ( z ) = ∑ k = 0 M − 1 z − k ( ∑ m = 0 K − 1 h ( k + m M ) z − m M ) = ∑ k = 0 M − 1 z − k H k ( z M ) H(z) = \sum_{k=0}^{M-1} z^{-k} \left( \sum_{m=0}^{K-1} h(k + mM) z^{-mM} \right) = \sum_{k=0}^{M-1} z^{-k} H_k(z^M) H(z)=k=0∑M−1​z−k(m=0∑K−1​h(k+mM)z−mM)=k=0∑M−1​z−kHk​(zM)

3 时域操作等价性

原FIR输出： y ( n ) = ∑ m = 0 N − 1 h ( m ) x ( n − m ) y(n) = \sum_{m=0}^{N-1} h(m)x(n-m) y(n)=m=0∑N−1​h(m)x(n−m) 多相结构输出： y ( n ) = ∑ k = 0 M − 1 ∑ m = 0 K − 1 h k ( m ) x ( n − k − m M ) y(n) = \sum_{k=0}^{M-1} \sum_{m=0}^{K-1} h_k(m) x\left(n - k - mM\right) y(n)=k=0∑M−1​m=0∑K−1​hk​(m)x(n−k−mM)

4、物理意义验证

1. 时域解释 原FIR的输出 y ( n ) = ∑ m = 0 N − 1 h ( m ) x ( n − m ) y(n) = \sum_{m=0}^{N-1} h(m)x(n-m) y(n)=∑m=0N−1​h(m)x(n−m)，等效于：

• 将输入 x ( n ) x(n) x(n)分为 M M M个相位分量（ x ( n − k ) x(n-k) x(n−k)， k = 0 , 1 , … , M − 1 k=0,1,\dots,M-1 k=0,1,…,M−1）
• 每个子滤波器 H k H_k Hk​处理对应的相位分量
• 结果相加得到输出 y ( n ) y(n) y(n)

2. 频域验证 通过替换 z = e j ω z = e^{j\omega} z=ejω，可验证原频率响应与多相分解后的响应一致。

5、应用场景

1. 多速率信号处理 在抽取（Decimation）和插值（Interpolation）中，多相分解可降低计算复杂度。
2. 并行化实现 各子滤波器 H k ( z M ) H_k(z^M) Hk​(zM)可并行计算，提升硬件效率。
3. 滤波器组设计 用于均匀DFT滤波器组、小波变换等。

二、Python实现与验证

该实战，通过两种同的方法进行抽取滤波，将信号采样率降低到原来的1/4，并对结果进行对比和误差分析。

1 生成测试信号

import numpy as np import matplotlib.pyplot as plt from scipy.signal import firwin, lfilter

生成测试信号：10Hz正弦波 + 100Hz高频噪声

fs = 1000 # 采样率 t = np.arange(0, 1, 1/fs) x = np.sin(2np.pi10t) + 0.5np.random.randn(len(t)) # 原始信号

2 设计FIR滤波器

设计低通FIR滤波器（截止频率50Hz，阶数N=31）

N = 31 fc = 50 h = firwin(N, fc, fs=fs, window=‘hamming’) # 获取系数

3 标准FIR滤波

y_fir = lfilter(h, 1, x) # FIR滤波结果

4 多相分解（M=4）

M = 4 # 分解因子 poly_h = [h[k::M] for k in range(M)] # 分解为M个子滤波器

补零对齐长度（确保所有子滤波器长度一致）

max_len = max(len(p) for p in poly_h) poly_h = [np.pad(p, (0, max_len - len(p))) for p in poly_h]

5 多相滤波实现

初始化输出

y_poly = np.zeros_like(x)

处理每个子滤波器分支

for k in range(M):

抽取输入信号的相位分量

x_k = x[k::M]

子滤波器滤波

y_k = lfilter(poly_h[k], 1, x_k)

上采样并添加延迟

y_up = np.zeros(len(x)) y_up[k::M] = y_k # 上采样（插入零） y_up = np.roll(y_up, -k) # 延迟补偿 y_poly += y_up # 合并分支结果

截断初始延迟

y_poly = y_poly[:len(x)-N] y_fir = y_fir[:len(x)-N] x_trim = x[:len(x)-N] y_fir = np.roll(y_fir, -3)

6 可视化对比

plt.figure(figsize=(12, 8))

原始信号与滤波结果

plt.subplot(3,1,1) plt.plot(x_trim, label=‘原始信号’, alpha=0.5) plt.plot(y_fir, label=‘FIR滤波结果’, linewidth=2) plt.legend() plt.title(‘FIR滤波效果’) y_fir = y_fir[0::M] #抽取 y_poly = y_poly[0::M] #抽取

多相滤波结果对比

plt.subplot(3,1,2) plt.plot(y_poly, label=‘多相滤波结果’, linestyle=’–’) plt.plot(y_fir, label=‘FIR滤波结果’, alpha=0.7) plt.legend() plt.title(‘多相滤波与FIR结果对比’)

误差分析

plt.subplot(3,1,3) error = y_fir - y_poly[:len(y_fir)] plt.plot(error, label=‘误差’, color=‘red’) plt.legend() plt.title(‘误差曲线 (最大误差: {:.2e})’.format(np.max(np.abs(error)))) plt.tight_layout() plt.show()

三、代码输出验证

1. 图形对比

• 第一张图展示原始信号与FIR滤波结果。
• 第二张图叠加显示FIR与多相滤波结果，两者应完全重合。
• 第三张图显示误差曲线。

2. 数值验证 误差曲线最大值接近机器精度，证明两种结构数学等价。

四、关键点说明

1. 多相分解的物理意义

• 每个子滤波器处理信号的特定相位分量，通过并行化降低计算复杂度。

2. 延迟补偿的重要性

• 分支信号需通过np.roll对齐时间轴，确保相位同步。

3. 应用场景优势

• 在多速率系统中（如抽取/插值），多相分解可减少计算量达 M M M倍。

通过上述实现，可直观验证FIR滤波器与其多相分解形式的等效性，为信号处理系统优化提供可靠依据。