📚 课程目录第2章基于互相关的时延估计

2. 基于互相关的时延估计

任务目标

本案例实现基于互相关的时延估计：给定两路含噪声的音频信号（其中一路相对另一路存在已知延迟），通过计算互相关函数并搜索峰值位置来估计延迟量。同时对比矩形窗与汉宁窗对估计精度的影响。

环境依赖

Python 3.8+，仅需NumPy和Matplotlib（用于可视化验证）。无需额外音频文件，使用合成信号。

核心实现

import numpy as np

def generate_signals(true_delay, n_samples=4096, fs=8000, noise_std=0.05):
    """生成两路含噪信号，x2相对x1延迟true_delay个采样点"""
    t = np.arange(n_samples) / fs
    freq = 440.0  # A4音符
    x1 = np.sin(2 * np.pi * freq * t)
    # 延迟通过补零实现（非整数延迟用插值，此处仅整数延迟）
    x2 = np.zeros_like(x1)
    if true_delay > 0:
        x2[true_delay:] = x1[:-true_delay]
    else:
        x2[:true_delay] = x1[-true_delay:]
    # 加入独立高斯白噪声
    rng = np.random.default_rng(42)
    x1 += rng.normal(0, noise_std, n_samples)
    x2 += rng.normal(0, noise_std, n_samples)
    return x1, x2

def cross_correlation(x1, x2, window=None):
    """计算互相关，返回延迟索引数组与对应相关值"""
    n = len(x1)
    if window is not None:
        w = window(n)
        x1 = x1 * w
        x2 = x2 * w
    # 使用FFT加速互相关计算
    fft_size = 2 * n
    X1 = np.fft.rfft(x1, fft_size)
    X2 = np.fft.rfft(x2, fft_size)
    # 互相关 = IFFT(X1 * conj(X2))
    corr = np.fft.irfft(X1 * np.conj(X2), fft_size)
    # 只保留有效延迟范围 [-n/2, n/2)
    lags = np.arange(-n//2, n//2)
    corr = np.concatenate([corr[-n//2:], corr[:n//2]])
    return lags, corr

def estimate_delay(x1, x2, window_fn=None):
    """返回估计延迟（采样点）"""
    lags, corr = cross_correlation(x1, x2, window_fn)
    peak_idx = np.argmax(np.abs(corr))
    return lags[peak_idx]

def evaluate_windows(true_delay=50, trials=100):
    """比较矩形窗与汉宁窗在不同噪声下的估计精度"""
    results = {'rect': [], 'hann': []}
    for _ in range(trials):
        x1, x2 = generate_signals(true_delay)
        d_rect = estimate_delay(x1, x2, window_fn=None)  # 矩形窗=不加窗
        d_hann = estimate_delay(x1, x2, window_fn=np.hanning)
        results['rect'].append(d_rect)
        results['hann'].append(d_hann)
    return results

# 单次演示
x1, x2 = generate_signals(true_delay=50)
d_rect = estimate_delay(x1, x2)
d_hann = estimate_delay(x1, x2, window_fn=np.hanning)
print(f"真实延迟: 50 采样点")
print(f"矩形窗估计: {d_rect} 采样点")
print(f"汉宁窗估计: {d_hann} 采样点")

# 统计评估（100次试验）
res = evaluate_windows()
rect_err = np.abs(np.array(res['rect']) - 50)
hann_err = np.abs(np.array(res['hann']) - 50)
print(f"\n矩形窗: 平均误差 = {rect_err.mean():.2f}, 正确率 = {(rect_err==0).mean()*100:.1f}%")
print(f"汉宁窗: 平均误差 = {hann_err.mean():.2f}, 正确率 = {(hann_err==0).mean()*100:.1f}%")

运行与验证

直接运行上述脚本，观察单次估计输出与统计评估结果。
判定标准：矩形窗在低噪声下应能正确估计延迟（误差为0）；汉宁窗因抑制了信号边缘不连续，在相同噪声下通常具有更高的正确估计率或更低的平均误差。
若需可视化互相关波形，可添加：import matplotlib.pyplot as plt; lags, corr = cross_correlation(x1, x2); plt.plot(lags, corr); plt.xlabel('延迟(采样点)'); plt.ylabel('互相关值'); plt.show()

常见问题

延迟为小数时：本实现仅支持整数采样点延迟。若需亚采样精度，需对互相关峰值附近进行抛物线插值或使用相位法。
窗函数的作用：矩形窗等效于直接截断，信号首尾不连续会产生频谱泄漏；汉宁窗将两端平滑衰减至零，减少泄漏但主瓣变宽，可能降低对邻近峰值的分辨能力。
噪声影响：当噪声标准差超过0.2时，两种窗的估计误差都会显著增大，此时可增加信号长度或使用加权互相关（如PHAT）改善。