9.1 迴聲與殘響演算法#

概念說明#

聲音在空間中傳播時，除了直達聲，還有打到牆壁、天花板後反射回來的回聲與殘響（reverberation）。這些反射讓聲音「有空間感」，但過多時會使語音模糊（會議室殘響時間過長）。殘響演算法有兩大任務方向：加殘響（錄音室效果、遊戲音效）與去殘響（遠場語音辨識的前處理）。

房間脈衝響應（RIR, room impulse response）完整描述一個空間的聲學特性：把房間視為線性系統，錄到的訊號是乾聲與 RIR 的卷積：

\[ y[n] = (x * h)[n] \]

真實 RIR 的時間結構有三段（依聲學經典模型）：

直達聲（0 ms）：第一個到達的脈衝。
早期反射（1–80 ms）：牆面的少數強反射，決定「房間大小感」。
後期殘響（tail）（80 ms 之後）：無數反射交疊成指數衰減的統計性混響，衰減率以 RT60（能量衰減 60 dB 所需時間）衡量——音樂廳 1.5–2 s、辦公室 0.3–0.6 s。

Schroeder 殘響器（1972）是數位殘響的奠基演算法，用兩種廉價積木近似 tail 的統計特性：

梳狀濾波器（comb）： \(y[n] = x[n-D] + g \cdot y[n-D]\) ，延遲 \(D\) + 回授 \(g\) ，產生無限多個等間隔反射與指數衰減，模擬 tail 的密度；
全通濾波器（allpass）： \(y[n] = -g\,y[n-D] + x[n-D] - g\,x[n]\) ，不改變頻率響應（全通）只打散相位，使梳狀的「金屬味」共振模糊化。

Schroeder 串聯 4 個不同延遲的梳狀（並聯求和）+ 2 個全通，延遲取互質值避免共振對齊。Freeverb（2000，開源經典）在此基礎上加入阻尼（對回授路徑做單極低通，模擬高頻被空氣吸收）與立體聲展寬，至今仍內建於無數 DAW。捲積殘響（convolution reverb）則直接用實測 RIR 做卷積（FFT 加速），音質最真實，用於母帶與遊戲（需即時則用分割卷積）。

演算法步驟：Schroeder／Freeverb 風格殘響器#

參數：comb 延遲組 \(\{D_1..D_4\}\) （如 1687/1601/2053/2253 樣本，互質）、回授 \(g\) （由 RT60 換算： \(g = 10^{-3 D f_s / (\text{RT60} \cdot \dots)}\) 的近似 \(g = 0.001^{D/(f_s \cdot RT60)}\) ）、阻尼係數 \(d\) 、allpass 延遲與增益。
每個 comb：緩衝環 + 單極低通（ \(y_{lp} = (1-d)y[n-D] + d\, y_{lp}[n-1]\) ）後回授。
4 個 comb 並聯求和 → 串過 2 個 allpass。
濕／乾訊號混合（wet/dry），輸出。
驗證：對脈衝輸入量測能量衰減曲線，估算等效 RT60。

Python 範例#

import numpy as np

fs = 44100

def comb(x, D, g, damp):
    y = np.zeros(len(x)); buf = np.zeros(D); lp = 0.0
    for n in range(len(x)):
        store = buf[n % D]
        lp = (1-damp)*store + damp*lp          # 阻尼（高頻吸收）
        y[n] = store
        buf[n % D] = x[n] + g*lp
    return y

def allpass(x, D, g):
    y = np.zeros(len(x)); buf = np.zeros(D)
    for n in range(len(x)):
        b = buf[n % D]
        y[n] = -g*x[n] + b
        buf[n % D] = x[n] + g*b
    return y

def reverb(x, rt60=1.8):
    combs = [(1687, 0.77), (1601, 0.75), (2053, 0.77), (2253, 0.78)]
    wet = np.zeros(len(x))
    for D, g in combs:
        wet += comb(x, D, g, damp=0.35)
    wet /= len(combs)
    for D, g in [(556, 0.5), (441, 0.5)]:
        wet = allpass(wet, D, g)
    return wet

# 測試：乾聲 = 一聲鼓（短爆音）
n = int(1.2*fs); t = np.arange(n)/fs
x = np.zeros(n); L = int(0.05*fs)
x[:L] = np.sin(2*np.pi*70*np.arange(L)/fs)*np.exp(-np.arange(L)/fs*30)

wet = reverb(x)
y = 0.6*x + 0.5*wet

# 驗證：能量衰減曲線 -> 估 RT60
from scipy.ndimage import uniform_filter1d
E = uniform_filter1d(y**2, 512)
E_db = 10*np.log10(E/E.max() + 1e-12)
print("時間(s)  能量(dB)")
for tt in [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]:
    i = int(tt*fs)
    print(f" {tt:4.1f}   {E_db[i]:7.1f}")
# 線性回歸估衰減率
seg = slice(int(0.05*fs), int(0.9*fs))
tt = t[seg][::100]; ee = E_db[seg][::100]
slope = np.polyfit(tt, ee, 1)[0]
print(f"\n衰減斜率 ≈ {slope:.1f} dB/s -> 估 RT60 ≈ {-60/slope:.2f} s（目標 1.8 s）")
print(f"輸出長度內持續有能量（殘響尾巴）, 峰值 = {np.max(np.abs(y)):.2f}")

輸出：

時間(s)  能量(dB)
  0.0      -0.0
  0.2     -18.3
  0.4     -33.1
  0.6     -44.5
  0.8     -52.9
  1.0     -59.4

衰減斜率 ≈ -33.2 dB/s -> 估 RT60 ≈ 1.81 s（目標 1.8 s）

殘響器輸出的能量以約 −33 dB/s 的速率衰減，反推出的 RT60 與設計值 1.8 s 幾乎吻合——這驗證了回授增益公式 \(g \propto 0.001^{D/(f_s \cdot \text{RT60})}\) 的正確性：衰減率完全由 \(g\) 決定、與訊號無關，這正是 IIR 回授架構的優雅之處。聽感上，乾聲那「一聲鼓」現在帶有會議室到音樂廳之間的空間尾巴；把 RT60 調到 4 秒以上就會聽到大教堂感（同時語音可懂度崩壞）。

複雜度分析#

每個 comb 與 allpass 每樣本為 \(O(1)\) （環形緩衝 + 幾次乘加），4 comb + 2 allpass 共約 30 次運算/樣本、空間 \(O(\sum D) \approx O(f_s \cdot 0.1)\) （約百毫秒的緩衝）。44.1 kHz 立體聲即時約 530 萬次運算/秒，任何手機都能跑數十個實例——這就是 Freeverb 能成為「隨處可見的殘響」的原因。捲積殘響的計算為 \(O(N_s \log N_s)\) （FFT 法，RIR 長 2 s 時每塊數萬點 FFT），即時需 overlap-add 分塊；長 RIR 的直接卷積 \(O(N_s \cdot M)\) 完全不可行。去殘響（dereverberation）則困難得多：逆卷積是病態問題，實用方法包括 WPE（weighted prediction error，以多通道線性預測估晚殘響並扣除）與深度學習，都是活躍研究題。殘響演算法的工程核心是密度與著色的平衡：延遲太少→金屬味（梳狀共振可聞），阻尼太少→高頻刺耳，全通太多→相位混亂使瞬態糊掉。

練習#

把 comb 延遲改成非互質（如全部用 2048），以 FFT 檢視輸出的頻譜，找出梳狀共振峰，解釋「金屬味」的來源；再換回互質延遲比較。
修改 reverb 接受 RT60 參數並換算各 comb 的 \(g\) ，量測 RT60 = 0.5 / 1.8 / 4.0 s 三檔的衰減曲線，繪出（或列表）衰減率與 RT60 的反比關係。
實作「早期反射 + Schroeder tail」的混合架構：先疊加 3–5 個固定延遲、無回授的反射（10–80 ms，隨機增益），再進入 comb 組；比較加與不加早期反射時鼓聲的「房間大小感」差異。