---

昨天的《信號處理之插值、抽取與多項濾波》，已經介紹了插值抽取的多項濾率，今天詳細介紹多項濾波的數學推導，并附上實戰仿真代碼。

一、數學變換推導

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=?\frac{N}{M}?$（向上取整）。

將原系數按$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濾波器與其多相分解形式的等效性，為信號處理系統優化提供可靠依據。