2. 晚期混响功率谱密度估计

任务说明

本案例实现基于噪声估计框架的晚期混响功率谱密度（PSD）估计。输入为单通道含噪语音信号（含混响与加性噪声），输出为逐帧、逐频点的晚期混响PSD估计值，供后续维纳滤波或谱减法使用。核心思路是将晚期混响视为一种“有色噪声”，采用最小值跟踪（Minimum Statistics）方法估计其时变PSD，并输出可直接用于抑制模块的估计结果。

算法原理与实现

采用Martin提出的最小值统计噪声估计框架：对带噪语音功率谱进行平滑，在固定长度窗口内跟踪平滑功率谱的最小值，用偏差补偿因子修正最小值得到噪声（此处即晚期混响）PSD估计。算法假设语音活动期间功率谱会显著高于纯混响段，因此最小值能够反映混响基底。实现中采用递归平滑和分块最小值搜索，兼顾跟踪速度与稳定性。

import numpy as np

def estimate_late_reverb_psd(x, fs=16000, nfft=512, hop=128,
                             alpha_s=0.7, win_len=125, alpha_p=0.2):
    """
    基于最小值统计的晚期混响PSD估计
    参数:
        x: 单通道时域信号 (含混响+噪声)
        fs: 采样率
        nfft: FFT长度
        hop: 帧移
        alpha_s: 平滑因子 (0~1)
        win_len: 最小值搜索窗口帧数
        alpha_p: 偏差补偿因子
    返回:
        PSD_rev: 晚期混响PSD估计 (帧数 x nfft/2+1)
        frames: 分帧后的复数谱 (帧数 x nfft/2+1)
    """
    # 分帧加窗
    frame_len = nfft
    n_frames = (len(x) - frame_len) // hop + 1
    frames = np.zeros((n_frames, nfft), dtype=np.complex128)
    win = np.hanning(frame_len)
    for i in range(n_frames):
        start = i * hop
        seg = x[start:start+frame_len] * win
        frames[i] = np.fft.fft(seg, nfft)

    # 功率谱
    spec = np.abs(frames[:, :nfft//2+1])**2
    n_bins = spec.shape[1]

    # 递归平滑
    P_smooth = np.zeros_like(spec)
    P_smooth[0] = spec[0]
    for t in range(1, n_frames):
        P_smooth[t] = alpha_s * P_smooth[t-1] + (1 - alpha_s) * spec[t]

    # 最小值跟踪 (分块)
    P_min = np.zeros_like(spec)
    n_blocks = n_frames // win_len
    for b in range(n_blocks):
        start = b * win_len
        end = min((b+1)*win_len, n_frames)
        block_min = np.min(P_smooth[start:end], axis=0)
        P_min[start:end] = block_min

    # 处理剩余帧
    if n_blocks * win_len < n_frames:
        start = n_blocks * win_len
        block_min = np.min(P_smooth[start:], axis=0)
        P_min[start:] = block_min

    # 偏差补偿
    B_min = 1.0 / (1.0 + alpha_p * (n_bins - 1))
    PSD_rev = B_min * P_min

    return PSD_rev, frames

# 示例：生成模拟输入 (混响语音+噪声)
np.random.seed(42)
fs = 16000
t = np.arange(fs * 2) / fs
clean = 0.5 * np.sin(2 * np.pi * 300 * t) + 0.3 * np.sin(2 * np.pi * 800 * t)
reverb = np.convolve(clean, np.exp(-np.arange(2000)/1000))[:len(clean)]  # 模拟混响
noise = 0.02 * np.random.randn(len(clean))
x = reverb + noise

# 运行估计
PSD_rev, frames = estimate_late_reverb_psd(x, fs=fs)

# 输出关键参数供验证
print(f"输入长度: {len(x)} 样本, 帧数: {PSD_rev.shape[0]}, 频点数: {PSD_rev.shape[1]}")
print(f"晚期混响PSD范围: {PSD_rev.min():.2e} ~ {PSD_rev.max():.2e}")
print(f"首帧前10个频点PSD: {PSD_rev[0, :10]}")

验证方法

运行上述代码，观察输出：

• 程序正常执行无报错，输出帧数与频点数符合预期（帧数≈(len(x)-nfft)//hop+1，频点=nfft//2+1=257）；
• PSD估计值均为非负有限数，且数值范围合理（与输入信号功率量级一致）；
• 在纯噪声段（如信号末尾0.5秒无语音激励部分），PSD估计应接近该段实际噪声功率；
• 可通过将PSD_rev与输入带噪语音功率谱对比，确认估计值低于语音活动段的总功率，高于静音段底噪，符合“晚期混响基底”的物理意义。

判定标准：代码运行无异常，输出维度正确，PSD值非负且处于合理量级，即视为通过。

落地与使用

本案例代码为独立Python脚本，依赖numpy。将代码保存为reverb_psd_est.py，在Python 3.8+环境执行python reverb_psd_est.py即可运行。实际应用时，将x替换为麦克风采集的含混响语音，调整win_len（建议100~200帧，对应1~2秒）以适应混响时间。估计得到的PSD_rev可直接作为维纳滤波器中噪声PSD项，或用于谱减法的过减因子控制。