信号处理抽取多项滤波的数学推导与仿真
昨天的《信号处理之插值、抽取与多项滤波》，已经介绍了插值抽取的多项滤率，今天详细介绍多项滤波的数学推导，并附上实战仿真代码。
一、数学变换推导
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−1h(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−1z−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−1h(n)z−n=k=0∑M−1m=0∑K−1h(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−1hk(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−1z−k(m=0∑K−1h(k+mM)z−mM)=k=0∑M−1z−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−1h(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−1m=0∑K−1hk(m)x(n−k−mM)
4、物理意义验证
-
时域解释
原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−1h(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)
-
频域验证
通过替换 z = e j ω z = e^{j\omega} z=ejω，可验证原频率响应与多相分解后的响应一致。
5、应用场景
- 多速率信号处理
在抽取（Decimation）和插值（Interpolation）中，多相分解可降低计算复杂度。 - 并行化实现
各子滤波器 H k ( z M ) H_k(z^M) Hk(zM)可并行计算，提升硬件效率。 - 滤波器组设计
用于均匀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(2*np.pi*10*t) + 0.5*np.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()
三、代码输出验证
-
图形对比
- 第一张图展示原始信号与FIR滤波结果。
- 第二张图叠加显示FIR与多相滤波结果，两者应完全重合。
- 第三张图显示误差曲线。
-
数值验证
误差曲线最大值接近机器精度，证明两种结构数学等价。
四、关键点说明
-
多相分解的物理意义
- 每个子滤波器处理信号的特定相位分量，通过并行化降低计算复杂度。
-
延迟补偿的重要性
- 分支信号需通过
np.roll对齐时间轴，确保相位同步。
- 分支信号需通过
-
应用场景优势
- 在多速率系统中（如抽取/插值），多相分解可减少计算量达 M M M倍。
通过上述实现，可直观验证FIR滤波器与其多相分解形式的等效性，为信号处理系统优化提供可靠依据。