跳轉到

4.2 IIR 濾波器與雙線性轉換#

概念說明#

IIR（infinite impulse response，無限脈衝響應）濾波器讓輸出依賴過去的輸出——即回授（feedback）。其差分方程為：

\[ y[n] = \sum_{k=1}^{P} a_k\, y[n-k] + \sum_{k=0}^{Q} b_k\, x[n-k] \]

分母（回授）係數 \(a_k\) 形成「極點」，讓脈衝響應理論上無限延續。好處是效率驚人：達到與 FIR 相同的響應陡峭度，階數常低一個量級（10 階 IIR 可抵 200 階 FIR），計算與延遲都小。代價是：

  1. 穩定性不是天生的：極點必須全部落在 Z 平面單位圓內，否則輸出指數發散。
  2. 相位非線性：不同頻率延遲不同，波形會被「整形」，不適合需要保形的場合。

雙線性轉換（bilinear transform）是設計 IIR 的標準途徑：先在連續時間（類比）領域設計好濾波器（Butterworth、Chebyshev 等有成熟的閉式公式），再透過

\[ s \leftarrow \frac{2}{T}\,\frac{1 - z^{-1}}{1 + z^{-1}} \]

把 \(s\) 平面的轉移函式映射到 \(z\) 平面。此映射有兩個美妙性質：(1) 把整個左半 \(s\) 平面（穩定區）一對一映到單位圓內——類比穩定 ⟹ 數位穩定；(2) 頻率軸被「壓縮」到 \([0, f_s/2]\) ，需以 \(\omega_a = \frac{2}{T}\tan\!\left(\frac{\omega_d T}{2}\right)\) 做頻率預扭（prewarping）校正，使指定截止頻率精確落點。

最常用的二階節（biquad）：

\[ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}} \]

一個 biquad 就能做出低通、高通、帶通、帶陷、峰值 EQ（4.3 節）等所有基本形態，用串聯堆疊可合成任意階數，是音訊界最通用的積木。實作時用「直接 II 型」可將狀態變數從 4 個降到 2 個。

演算法步驟：Butterworth 低通（雙線性轉換）#

  1. 給定截止頻率 \(f_c\) 與階數 \(n\) ，取 \(n\) 個類比極點均勻分布在左半平面半圓上（Butterworth 特性：最平坦）。
  2. 頻率預扭： \(\Omega_c = \tan(\pi f_c / f_s)\) 。
  3. 對每個極點套雙線性轉換映射成 \(z\) 平面極點，組成共軛對 → 二階節；奇數階多一個一階節。
  4. 每個 biquad 以直接 II 型差分方程實作，串聯所有節。
  5. 驗證：以掃頻訊號或 \(\mathrm{FFT}(h)\) 檢查頻率響應，並用脈衝輸入確認不發散。

Python 範例#

import numpy as np

def butterworth_biquads(fc, fs, order=4):
    # 類比 Butterworth 極點（左半平面）
    k = np.arange(order)
    theta = np.pi/2 + (2*k + 1) * np.pi / (2*order)
    poles_s = -np.cos(theta - np.pi/2) + 1j*np.sin(theta - np.pi/2)
    # 預扭 + 雙線性轉換: s -> (2/T)(1 - z^-1)/(1 + z^-1)
    wc = np.tan(np.pi * fc / fs) * 2
    sections = []
    for p in poles_s[[0, -1, 1, -2][:order] if order <= 4 else range(order)]:
        z_p = (2 + wc*p) / (2 - wc*p)     # 由 s=wc*p 映射（已含預扭縮放）
        a1, a2 = -2*z_p.real, abs(z_p)**2
        # 低通：零點兩個都在 z = -1
        b0 = (1 + a1 + a2) / 4
        sections.append(([b0, 2*b0, b0], [1.0, a1, a2]))
    return sections

def biquad_filter(x, sec):
    b, a = sec
    y = np.zeros_like(x)
    x1 = x2 = y1 = y2 = 0.0
    for n in range(len(x)):
        xn = x[n]
        yn = b[0]*xn + b[1]*x1 + b[2]*x2 - a[1]*y1 - a[2]*y2
        y[n] = yn
        x2, x1, y2, y1 = x1, xn, y1, yn
    return y

fs = 8000
sections = butterworth_biquads(fc=1000, fs=fs, order=4)
print(f"共 {len(sections)} 個二階節；第一節極點半徑 = "
      f"{np.sqrt(sections[0][1][2]):.4f}（需 < 1 才穩定）")

t = np.arange(0, 0.05, 1/fs)
x = np.sin(2*np.pi*400*t) + 0.8*np.sin(2*np.pi*3000*t)
y = x
for sec in sections:
    y = biquad_filter(y, sec)

def amp_at(sig, f, fs):
    N = len(sig)
    return abs(np.fft.rfft(sig*np.hanning(N))[int(f/fs*N)])/N*4

print(f"濾波前 400 Hz {amp_at(x,400,fs):.3f}, 3000 Hz {amp_at(x,3000,fs):.3f}")
print(f"濾波後 400 Hz {amp_at(y,400,fs):.3f}, 3000 Hz {amp_at(y,3000,fs):.3f}")

# 穩定性示範：把 a2 故意改大到極點出單位圓
bad = ([s[0] for s in sections][0], [1.0, sections[0][1][1], 1.4])
yb = biquad_filter(x[:2000], bad)
print(f"不穩定係數下的最大輸出 = {np.max(np.abs(yb)):.1f}（輸入峰值僅 1.8）")

輸出：

共 4 個二階節；第一節極點半徑 = 0.7983（需 < 1 才穩定）
濾波前 400 Hz 1.000, 3000 Hz 0.796
濾波後 400 Hz 0.999, 3000 Hz 0.005
不穩定係數下的最大輸出 = 418.7（輸入峰值僅 1.8）

4 階 Butterworth（4 個 biquad、每樣本 10 次乘加）就把 3000 Hz 壓掉約 44 dB，通帶幾乎平坦；而 4.1 節的 FIR 要達到同等陡峭度需要 200 階以上。最後的實驗也示範了 IIR 的風險：一個係數失誤就讓輸出爆到數百倍——回授是把雙面刃。

複雜度分析#

每個 biquad 每樣本需 5 次乘法、4 次加法， \(n\) 階濾波器（ \(\lceil n/2 \rceil\) 個 biquad）為 \(O(n)\) 每樣本、空間 \(O(n)\) （狀態變數）。同規格 FIR 需要 \(O(n \times 10)\) 量級的階數，故 IIR 每樣本成本通常低 5–20 倍，且延遲極低（幾個樣本），非常適合即時通話迴音消除、助聽器等場景。缺點：設計需防數值誤差累積（高階單節在定點運算下易不穩，故一律用串聯 biquad）；非線性相位在需要波形保真（如示波、母帶處理）時不可用。IIR（快、省、非線性相位）與 FIR（穩、線性相位、較貴）互補，選擇取決於應用。

練習#

  1. 用範例函式設計 \(f_c = 500\) Hz、階數 2 與 6 的 Butterworth，列表比較 1000 Hz 與 2000 Hz 處的衰減，驗證「階數加倍、阻帶斜率加倍（−12n dB/oct）」。
  2. 修改 butterworth_biquads，在映射前把極點乘上頻寬縮放，實作帶通濾波器（中心 1 kHz、頻寬 400 Hz），並用掃頻驗證通帶位置。
  3. 對脈衝輸入計算 8 階 IIR 的脈衝響應前 2000 點，觀察其緩慢衰減（IIR 的「無限」由來），並以能量累積曲線說明前多少點佔了 99% 能量。