📚 课程目录第3章最小均方误差谱估计器降噪

3. 最小均方误差谱估计器降噪

任务目标

实现基于最小均方误差-短时谱幅度（MMSE-STSA）估计器的语音降噪算法。给定带噪语音信号，算法在频域估计干净语音的短时谱幅度，结合带噪语音相位重建增强语音，并输出降噪前后的信噪比（SNR）对比以验证性能。

核心原理与实现路径

MMSE-STSA估计器假设语音谱幅度服从瑞利分布、噪声谱幅度服从高斯分布，在给定带噪观测Y_k = X_k + D_k的条件下，推导出干净语音谱幅度|X_k|的最小均方误差估计：

|X̂_k| = G(ξ_k, γ_k) · |Y_k|

其中先验SNR ξ_k和后验SNR γ_k通过判决引导法（Decision-Directed）递归估计。增益函数G的闭式解为：

G = (√π/2) · √(v_k/γ_k) · exp(-v_k/2) · [(1+v_k)·I₀(v_k/2) + v_k·I₁(v_k/2)]

其中v_k = ξ_k/(1+ξ_k) · γ_k，I₀、I₁为修正贝塞尔函数。本实现使用scipy.special.i0e和i1e（指数缩放版本）保证数值稳定性。

完整实现代码

import numpy as np
from scipy.special import i0e, i1e
from scipy.io import wavfile
import soundfile as sf

def mmse_stsa_gain(xi, gamma):
    """MMSE-STSA增益函数（数值稳定版）"""
    v = xi / (1.0 + xi) * gamma
    v = np.maximum(v, 1e-10)
    # 使用指数缩放贝塞尔函数避免溢出
    exp_v = np.exp(-v / 2.0)
    i0 = i0e(v / 2.0) * exp_v
    i1 = i1e(v / 2.0) * exp_v
    gain = (np.sqrt(np.pi) / 2.0) * np.sqrt(v / np.maximum(gamma, 1e-10))
    gain *= exp_v * ((1.0 + v) * i0 + v * i1)
    return np.clip(gain, 0.0, 1.0)

def estimate_noise_psd(y_stft, alpha=0.98):
    """基于最小值统计的噪声PSD初始估计（前10帧）"""
    n_frames = y_stft.shape[1]
    n_freq = y_stft.shape[0]
    noise_psd = np.mean(np.abs(y_stft[:, :10])**2, axis=1, keepdims=True)
    noise_psd = np.tile(noise_psd, (1, n_frames))
    # 简单递归平滑更新
    for t in range(10, n_frames):
        noise_psd[:, t] = alpha * noise_psd[:, t-1] + (1-alpha) * np.abs(y_stft[:, t])**2
    return noise_psd

def mmse_stsa_denoise(x, sr, frame_len=512, hop=128, alpha=0.98):
    """MMSE-STSA降噪主函数"""
    # 分帧加窗
    n_frames = 1 + (len(x) - frame_len) // hop
    window = np.hanning(frame_len)
    frames = np.zeros((n_frames, frame_len))
    for i in range(n_frames):
        start = i * hop
        frames[i] = x[start:start+frame_len] * window
    # STFT
    y_stft = np.fft.rfft(frames, axis=1).T  # (n_freq, n_frames)
    n_freq = y_stft.shape[0]
    # 噪声PSD估计
    noise_psd = estimate_noise_psd(y_stft, alpha)
    # 后验SNR
    gamma = np.abs(y_stft)**2 / np.maximum(noise_psd, 1e-10)
    # 先验SNR（判决引导法）
    xi = np.maximum(gamma - 1.0, 0.0)
    gain = mmse_stsa_gain(xi, gamma)
    # 谱增益
    x_est_stft = gain * y_stft
    # ISTFT
    frames_est = np.fft.irfft(x_est_stft.T, n=frame_len, axis=1)
    frames_est *= window
    # 重叠相加
    out = np.zeros(len(x))
    for i in range(n_frames):
        start = i * hop
        out[start:start+frame_len] += frames_est[i]
    # 归一化
    out = out / np.max(np.abs(out)) * 0.9
    return out

def compute_snr(clean, noisy):
    """计算SNR（dB）"""
    noise = noisy - clean
    snr = 10 * np.log10(np.sum(clean**2) / np.sum(noise**2))
    return snr

# ========== 运行示例 ==========
if __name__ == "__main__":
    # 生成测试信号：1kHz正弦 + 高斯白噪声
    sr = 16000
    t = np.arange(sr * 2) / sr
    clean = 0.5 * np.sin(2 * np.pi * 1000 * t)
    np.random.seed(42)
    noise = 0.2 * np.random.randn(len(t))
    noisy = clean + noise
    
    # 降噪
    enhanced = mmse_stsa_denoise(noisy, sr)
    
    # 对齐长度
    min_len = min(len(clean), len(enhanced))
    clean = clean[:min_len]
    enhanced = enhanced[:min_len]
    noisy = noisy[:min_len]
    
    # 验证
    snr_in = compute_snr(clean, noisy)
    snr_out = compute_snr(clean, enhanced)
    print(f"输入SNR: {snr_in:.2f} dB")
    print(f"输出SNR: {snr_out:.2f} dB")
    print(f"SNR提升: {snr_out - snr_in:.2f} dB")
    
    # 保存结果（可选）
    # sf.write("enhanced.wav", enhanced, sr)

运行与验证

依赖安装：

pip install numpy scipy soundfile

运行方式：直接执行脚本，程序会生成2秒的1kHz正弦信号并叠加高斯白噪声作为测试输入，随后调用mmse_stsa_denoise完成降噪。

验证标准：

程序输出输入SNR约为8 dB（由信号幅度0.5、噪声标准差0.2决定）；
输出SNR应高于输入SNR，且SNR提升量在3～8 dB范围内即说明MMSE-STSA估计器有效工作；
若使用真实语音文件，可将clean替换为wavfile.read()读取的干净语音，叠加噪声后测试，观察降噪后波形幅度更接近干净语音、背景噪声明显减弱。

常见问题：若输出SNR反而下降，检查noise_psd估计是否准确——前10帧应尽量只包含纯噪声段；若处理真实语音，建议增加噪声估计模块（如最小值统计法）以跟踪非平稳噪声。