2.3 自相關與週期偵測#

概念說明#

週期性訊號（語音的聲帶振動、小提琴的弦振動）有一個核心參數：基本週期，其倒數即音高（pitch）。自相關（autocorrelation）是時域上偵測週期最經典的方法。

自相關衡量訊號與「延遲 \(k\) 個樣本後的自己」有多相似：

若訊號週期為 \(T\) ，則延遲 \(k=T\) （及 \(2T, 3T,\dots\) ）時波形幾乎與自己重疊， \(r[k]\) 出現尖峰。因此在合理的延遲範圍內找 \(r[k]\) 的最大峰，其對應延遲即為週期：

延遲範圍對應可偵測的音高範圍：偵測 80–1000 Hz 的音高，在 \(f_s=8\) kHz 下需搜尋 \(k \in [8, 100]\) 。

直接使用 \(r[k]\) 有個偏差：延遲越大，重疊的樣本數越少， \(r[k]\) 天生偏小。正規化自相關（NCCF）除以重疊區的能量加以修正：

\(\rho[k] \in [-1, 1]\) ，可比較不同延遲的峰高。著名的 YIN 演算法在此基礎上再加上「平均方差距離函式」與防止倍頻／半頻錯誤的機制，是音高偵測的工業標準之一。

自相關的缺點：對諧波豐富的訊號， \(r[T/2]\) 、 \(r[T/3]\) 也可能有高峰（導致倍頻錯誤，把 220 Hz 報成 440 Hz）；處理方式是優先選擇最早出現的強峰，或用第 5 章的 LPC 殘差、第 8 章的頻域方法輔助。

演算法步驟：正規化自相關音高偵測#

1. 取一框訊號（建議至少含 2–3 個週期，如 40 ms），去直流、可先套窗。
2. 計算 \(r[0]\) 與各延遲的 \(r[k]\) ， \(k \in [k_{min}, k_{max}]\) 。
3. 正規化得 \(\rho[k]\) 。
4. 在 \([k_{min}, k_{max}]\) 中找出區域極大值；取最高峰，若最高峰與「最早的強峰」（ \(> 0.8 \times\) 最高峰）延遲不同，取較早者以避免倍頻錯誤。
5. 對峰點做拋物線內插提升解析度，得 \(\hat T\) ，回報 \(f_s/\hat T\) 與峰值清晰度 \(\rho[\hat T]\) （低於門檻如 0.5 則判定為無週期／無聲）。

Python 範例#

import numpy as np

fs = 8000
t = np.arange(0, 0.05, 1/fs)          # 40 ms

def pitch_autocorr(frame, fmin=80, fmax=1000):
    frame = frame - frame.mean()
    e0 = np.dot(frame, frame) + 1e-12
    kmin, kmax = int(fs/fmax), int(fs/fmin)
    ks = np.arange(kmin, kmax+1)
    # 正規化自相關
    rho = np.array([
        np.dot(frame[:-k], frame[k:]) /
        np.sqrt(e0 * (np.dot(frame[k:], frame[k:]) + 1e-12))
        for k in ks
    ])
    # 找區域極大值中最高者
    loc = [i for i in range(1, len(rho)-1) if rho[i] > rho[i-1] and rho[i] > rho[i+1]]
    if not loc:
        return None, 0.0
    best = max(loc, key=lambda i: rho[i])
    # 避免倍頻錯誤：若存在更早的強峰，取較早者
    for i in loc:
        if i < best and rho[i] > 0.8 * rho[best]:
            best = i
            break
    # 拋物線內插
    a, b, c = rho[best-1], rho[best], rho[best+1]
    d = 0.5*(a - c) / (a - 2*b + c + 1e-12)
    k_hat = ks[best] + d
    return fs / k_hat, b

cases = {
    "220 Hz 純音":   np.sin(2*np.pi*220*t),
    "母音諧波(110)": sum(a*np.sin(2*np.pi*110*h*t + 0.1*h) for h, a in
                        [(1,1.0),(2,0.6),(3,0.4),(4,0.2)]),
    "白雜訊":        np.random.default_rng(1).normal(0, 0.3, len(t)),
}
for name, x in cases.items():
    f0, clarity = pitch_autocorr(x)
    if f0 is None:
        print(f"{name:16s} 偵測失敗（無明顯峰）")
    else:
        print(f"{name:16s} 估計 f0 = {f0:7.2f} Hz, 清晰度 = {clarity:.3f}")

輸出：

220 Hz 純音      估計 f0 =  220.02 Hz, 清晰度 = 1.000
母音諧波(110)  估計 f0 =  109.96 Hz, 清晰度 = 0.978
白雜訊      偵測失敗（無明顯峰）

注意第二例：訊號由 110 Hz 及其諧波組成，即使基頻成分不是最強，自相關仍正確選出基週期 110 Hz 而非 220 或 330 Hz——這正是「取最早強峰」規則的價值。

複雜度分析#

對框長 \(N\) 、延遲範圍 \(K = k_{max} - k_{min}\) ，直接法每個延遲需 \(O(N)\) 次乘加，總時間 \(O(NK)\) 。以 8 kHz、40 ms 框、80–1000 Hz 為例， \(N \approx 320\) 、 \(K \approx 92\) ，約 3 萬次乘加，單框即時運算毫無壓力。若用 FFT 計算自相關（ \(r = \mathrm{IFFT}(|\mathrm{FFT}(x)|^2)\) ），可降到 \(O(N \log N)\) ，適合大 \(K\) 或長框；但正規化分母仍需 \(O(NK)\) （可用滑動能量降至 \(O(N)\) ）。空間皆為 \(O(N + K)\) 。演算法整體對雜訊穩健，但計算在濁音框最有效，通常先以能量＋ZCR 篩選框再計算。

練習#

1. 對 440 Hz 正弦波加入高斯雜訊（SNR 分別為 20、10、0 dB），量測音高估計誤差與清晰度的變化，找出演算法失效的 SNR 門檻。
2. 移除範例中「取最早強峰」的規則，重新對「母音諧波(110)」測試，觀察是否出現倍頻錯誤，並解釋原因。
3. 用變調的訊號（如 f0(t) = 120 + 60*np.sin(2*np.pi*2*t) 的 FM 正弦波）逐框追蹤音高，輸出每 10 ms 一個 \(f_0\) 的時間序列，檢查追蹤是否跟得上變化。