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

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

一、数学变换推导

1. 多相分解的核心思想

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

2. 数学变换推导

FIR滤波器的系统函数可表示为：
$$H(z)=\sum _{n=0}^{N-1}h(n)z^{-n}$$
其中，$h(n)$为滤波器系数，$N$为阶数。

设分解因子为$M$，则第$k$个子滤波器系数为：
$$h_k(m)=h(k+mM),0\le k<M$$

将FIR滤波器拆分为$M$个并行的子滤波器（即多相分量），每个子滤波器处理信号的特定相位分量，再通过延迟和组合实现等效。目标形式为：

$$H(z)=\sum _{k=0}^{M-1}{z}^{-k}{H}_k(z^M)$$
其中，$H_k(z^M)$表示第$k$个子滤波器的系统函数。

将$H(z)$按$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)}$$
其中，$K=\lceil \frac{N}{M}\rceil$（向上取整）。

将原系数按$M$个相位分组：

• 第$k$个子滤波器的系数为：$h_k(m)=h(k+mM)$
• 其系统函数为：
  $$H_k(z^M)=\sum _{m=0}^{K-1}{h}_k(m)z^{-mM}$$
  将$H(z)$重写为：
  $$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)$$

3. 时域操作等价性

原FIR输出：
$$y(n)=\sum _{m=0}^{N-1}h(m)x(n-m)$$
多相结构输出：
$$y(n)=\sum _{k=0}^{M-1}\sum _{m=0}^{K-1}{h}_k(m)x\left(n-k-mM\right)$$

4、物理意义验证

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

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

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

5、应用场景

1. 多速率信号处理
   在抽取（Decimation）和插值（Interpolation）中，多相分解可降低计算复杂度。
2. 并行化实现
   各子滤波器$H_k(z^M)$可并行计算，提升硬件效率。
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(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()

三、代码输出验证

1. 图形对比

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

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

四、关键点说明

1. 多相分解的物理意义

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

2. 延迟补偿的重要性

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

3. 应用场景优势

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

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