📚 课程目录第3章广义互相关-相位变换时延估计

3. 广义互相关-相位变换时延估计

任务目标

在混响环境下，基本互相关（CC）时延估计因多径效应导致峰值模糊。本案例实现广义互相关-相位变换（GCC-PHAT）算法，通过频域加权白化抑制混响影响，并与基本互相关方法对比估计精度。

核心原理

两麦克风接收信号模型：\(x_1[n]=s[n]*h_1[n]+v_1[n]\)，\(x_2[n]=s[n-D]*h_2[n]+v_2[n]\)。GCC-PHAT在频域计算互功率谱并归一化：

\[ \hat{R}_{12}(\tau)=\int_{-\infty}^{\infty} \frac{X_1(f)X_2^*(f)}{|X_1(f)X_2^*(f)|} e^{j2\pi f\tau} df \]

分母归一化将各频点幅度置为1，仅保留相位信息，从而抑制混响引起的幅度畸变。时延估计为 \(\hat{D}=\arg\max_\tau \hat{R}_{12}(\tau)\)。

环境依赖

需要：Python 3.8+，numpy，scipy（用于信号滤波与读取）。无第三方音频库依赖。

完整实现

import numpy as np
from scipy.signal import fftconvolve, lfilter

def generate_test_signals(fs=16000, T=0.5, D_true=12, rt60=0.3):
    """生成带混响的双通道测试信号（模拟房间冲激响应）"""
    n = int(fs * T)
    t = np.arange(n) / fs
    # 源信号：带通噪声激励
    rng = np.random.default_rng(42)
    s = lfilter([0.01], [1, -0.99], rng.standard_normal(n))
    # 模拟房间冲激响应（指数衰减+早期反射）
    h_len = int(rt60 * fs)
    h1 = np.exp(-np.arange(h_len) / (rt60 * fs / 6)) * rng.standard_normal(h_len)
    h2 = np.zeros(h_len + abs(D_true))
    h2[abs(D_true):] = np.exp(-np.arange(h_len) / (rt60 * fs / 6)) * rng.standard_normal(h_len)
    # 卷积生成接收信号
    x1 = fftconvolve(s, h1)[:n]
    x2 = fftconvolve(s, h2)[:n]
    # 添加传感器噪声
    x1 += 0.01 * rng.standard_normal(n)
    x2 += 0.01 * rng.standard_normal(n)
    return x1, x2, D_true

def gcc_phat(x1, x2, fs, max_lag=None):
    """GCC-PHAT时延估计，返回时延（秒）"""
    n = len(x1) + len(x2) - 1
    X1 = np.fft.rfft(x1, n)
    X2 = np.fft.rfft(x2, n)
    cross = X1 * np.conj(X2)
    # 相位变换加权
    cross_phat = cross / (np.abs(cross) + 1e-10)
    r = np.fft.irfft(cross_phat, n)
    # 搜索峰值范围
    if max_lag is None:
        max_lag = int(0.05 * fs)  # 默认搜索±50ms
    lags = np.arange(-max_lag, max_lag + 1)
    idx = (lags + n) % n
    peak_idx = lags[np.argmax(r[idx])]
    return peak_idx / fs

def basic_cc(x1, x2, fs, max_lag=None):
    """基本互相关时延估计"""
    n = len(x1) + len(x2) - 1
    X1 = np.fft.rfft(x1, n)
    X2 = np.fft.rfft(x2, n)
    cross = X1 * np.conj(X2)
    r = np.fft.irfft(cross, n)
    if max_lag is None:
        max_lag = int(0.05 * fs)
    lags = np.arange(-max_lag, max_lag + 1)
    idx = (lags + n) % n
    peak_idx = lags[np.argmax(r[idx])]
    return peak_idx / fs

# 主程序：生成信号并对比
fs = 16000
x1, x2, D_true = generate_test_signals(fs=fs, D_true=12)
delay_cc = basic_cc(x1, x2, fs)
delay_phat = gcc_phat(x1, x2, fs)

print(f"真实时延: {D_true/fs*1000:.2f} ms")
print(f"基本互相关估计: {delay_cc*1000:.2f} ms (误差 {abs(delay_cc-D_true/fs)*1000:.2f} ms)")
print(f"GCC-PHAT估计: {delay_phat*1000:.2f} ms (误差 {abs(delay_phat-D_true/fs)*1000:.2f} ms)")

运行与验证

将上述代码保存为 gcc_phat_demo.py，在Python环境中运行：python gcc_phat_demo.py
程序自动生成带混响的双通道信号（真实时延12采样点=0.75ms），分别用两种方法估计时延。
通过条件：GCC-PHAT估计误差应显著小于基本互相关误差（在rt60=0.3s混响下，GCC-PHAT误差通常<0.1ms，而基本互相关可能偏差数毫秒）。输出中GCC-PHAT误差值应小于CC误差值。

参数影响说明

rt60：混响时间（秒）。增大该值（如0.5~0.8）时，基本互相关性能急剧下降，而GCC-PHAT仍能保持较高精度，可自行修改验证。
D_true：真实时延（采样点）。修改后需同步调整输出对比。
若环境噪声过大（0.01系数调高），两种方法误差都会增大，但GCC-PHAT仍相对更稳健。

常见问题

峰值搜索范围：代码默认搜索±50ms，若实际时延超出此范围需调整max_lag参数。
频率分辨率：FFT长度取两信号长度之和减1，避免循环卷积混叠。
数值稳定性：分母添加1e-10防止除零。