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

---

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

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、物理意义验证

• 时域表明
  原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)
    
  
  
• 频域验证
  通过替换                                             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. 生成测试信号

1. import numpy as np
2. import matplotlib.pyplot as plt
3. from scipy.signal import firwin, lfilter
5. # 生成测试信号：10Hz正弦波 + 100Hz高频噪声
6. fs = 1000  # 采样率
7. t = np.arange(0, 1, 1/fs)
8. x = np.sin(2*np.pi*10*t) + 0.5*np.random.randn(len(t))  # 原始信号

复制代码

2. 设计FIR滤波器

1. # 设计低通FIR滤波器（截止频率50Hz，阶数N=31）
2. N = 31
3. fc = 50
4. h = firwin(N, fc, fs=fs, window='hamming')  # 获取系数

复制代码

3. 标准FIR滤波

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

复制代码

4. 多相分解（M=4）

1. M = 4  # 分解因子
2. poly_h = [h[k::M] for k in range(M)]  # 分解为M个子滤波器
4. # 补零对齐长度（确保所有子滤波器长度一致）
5. max_len = max(len(p) for p in poly_h)
6. poly_h = [np.pad(p, (0, max_len - len(p))) for p in poly_h]

复制代码

5. 多相滤波实现

1. # 初始化输出
2. y_poly = np.zeros_like(x)
4. # 处理每个子滤波器分支
5. for k in range(M):
6.     # 抽取输入信号的相位分量
7.     x_k = x[k::M]
8.    
9.     # 子滤波器滤波
10.     y_k = lfilter(poly_h[k], 1, x_k)
11.    
12.     # 上采样并添加延迟
13.     y_up = np.zeros(len(x))
14.     y_up[k::M] = y_k  # 上采样（插入零）
15.     y_up = np.roll(y_up, -k)  # 延迟补偿
16.    
17.     y_poly += y_up  # 合并分支结果
19. # 截断初始延迟
20. y_poly = y_poly[:len(x)-N]
21. y_fir = y_fir[:len(x)-N]
22. x_trim = x[:len(x)-N]
23. y_fir = np.roll(y_fir, -3)

复制代码

6. 可视化对比

1. plt.figure(figsize=(12, 8))
3. # 原始信号与滤波结果
4. plt.subplot(3,1,1)
5. plt.plot(x_trim, label='原始信号', alpha=0.5)
6. plt.plot(y_fir, label='FIR滤波结果', linewidth=2)
7. plt.legend()
8. plt.title('FIR滤波效果')
10. y_fir = y_fir[0::M] #抽取
11. y_poly = y_poly[0::M] #抽取
13. # 多相滤波结果对比
14. plt.subplot(3,1,2)
15. plt.plot(y_poly, label='多相滤波结果', linestyle='--')
16. plt.plot(y_fir, label='FIR滤波结果', alpha=0.7)
17. plt.legend()
18. plt.title('多相滤波与FIR结果对比')
20. # 误差分析
21. plt.subplot(3,1,3)
22. error = y_fir - y_poly[:len(y_fir)]
23. plt.plot(error, label='误差', color='red')
24. plt.legend()
25. plt.title('误差曲线 (最大误差: {:.2e})'.format(np.max(np.abs(error))))
27. plt.tight_layout()
28. plt.show()

复制代码

---

三、代码输出验证

• 图形对比
  
  
  
  • 第一张图展示原始信号与FIR滤波效果。
    
  • 第二张图叠加显示FIR与多相滤波效果，两者应完全重合。
    
  • 第三张图显示误差曲线。
    
  
  
• 数值验证
  误差曲线最大值接近机器精度，证明两种布局数学等价。
  

---

四、关键点阐明

• 多相分解的物理意义
  
  
  
  • 每个子滤波器处理信号的特定相位分量，通过并行化降低计算复杂度。
    
  
  
• 延迟赔偿的紧张性
  
  
  
  • 分支信号需通过np.roll对齐时间轴，确保相位同步。
    
  
  
• 应用场景上风
  
  
  
  • 在多速率体系中（如抽取/插值），多相分解可减少计算量达                                                  M                                          M                           M倍。
    
  
  

---

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

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。