谱减法（Spectral Subtraction）

谱减法是一种经典且广泛应用的语音增强算法，主要用于去除语音信号中的加性背景噪声。其核心思想非常简单：假设噪声是加性的且在语音活动间隙（只有噪声）的统计特性相对平稳，那么可以通过估计噪声的功率谱，并从带噪语音的功率谱中减去该噪声功率谱的估计值，从而得到增强后的语音功率谱。最后，结合带噪语音的相位信息（通常认为相位信息对听觉影响较小，故保留原相位），进行逆傅里叶变换得到增强后的时域语音信号。

核心步骤

1. 分帧与加窗： 将输入的带噪语音信号 $y(n)$ 分割成短时重叠帧，通常帧长为20-30ms，帧移为10ms。对每帧信号应用窗函数（如汉明窗）以减少频谱泄露。
2. 短时傅里叶变换（STFT）： 对每一帧信号进行STFT，得到其复数频谱 $Y(k,m)$，其中 $k$ 是频率索引，$m$ 是帧索引。
   $$Y(k,m)=\sum _{n=0}^{N-1}y(n+mH)w(n)e^{-j2\pi kn/N}$$
   这里 $N$ 是帧长，$H$ 是帧移，$w(n)$ 是窗函数。
3. 噪声功率谱估计： 在语音活动检测（VAD）判定为“静音”或“只有噪声”的帧上，估计噪声的平均功率谱 $|\hat{D}(k){|}^2$。最简单的方法是计算这些帧的幅度谱平方的平均值。
   $$|\hat{D}(k){|}^2=\frac{1}{M}\sum _{m\in \text{noise frames}}|Y(k,m){|}^2$$
   其中 $M$ 是用于估计的噪声帧数量。
4. 带噪语音功率谱估计： 计算当前帧带噪语音的功率谱 $|Y(k,m){|}^2$。
5. 谱减： 从带噪语音的功率谱中减去估计的噪声功率谱，得到增强语音的功率谱估计 $|\hat{X}(k,m){|}^2$。为了避免负值（这在物理上没有意义），通常采用半波整流或设置一个下限（如 $\beta |\hat{D}(k){|}^2$）。
   $$|\hat{X}(k,m){|}^2=\max \left(|Y(k,m){|}^2-\alpha |\hat{D}(k){|}^2,\beta |\hat{D}(k){|}^2\right)$$
   其中 $\alpha$ 是过减因子（通常 $\alpha \ge 1$，用于补偿噪声估计的不足）， $\beta$ 是谱下限因子（通常 $0<\beta \ll 1$，防止过度抑制）。
6. 增强幅度谱： 计算增强后的幅度谱 $|\hat{X}(k,m)|$。
   $$|\hat{X}(k,m)|=\sqrt{|\hat{X}(k,m){|}^2}$$
7. 相位保留与合成： 使用带噪语音的相位 $\angle Y(k,m)$ 和增强后的幅度谱 $|\hat{X}(k,m)|$ 重构复数频谱 $\hat{X}(k,m)$。
   $$\hat{X}(k,m)=|\hat{X}(k,m)|e^{j\angle Y(k,m)}$$
8. 逆STFT（ISTFT）： 对每一帧的 $\hat{X}(k,m)$ 进行逆STFT，得到增强后的时域信号帧 $\hat{x}(n,m)$。
   $$\hat{x}(n,m)=\frac{1}{N}\sum _{k=0}^{N-1}\hat{X}(k,m)e^{j2\pi kn/N}$$
9. 重叠相加（OLA）： 将处理后的帧通过重叠相加的方法合成最终的增强语音信号 $\hat{x}(n)$。

关键问题与改进

• 音乐噪声（Musical Noise）： 由于噪声估计误差和谱减的非线性操作，增强后的语音中可能出现类似音乐的随机尖峰（残留噪声），听起来像“啁啾”声或“嗡嗡”声。这是基本谱减法的主要缺点。
• 改进方法： 为了解决音乐噪声等问题，提出了多种改进方案：
  • 过减因子和谱下限： 合理选择 $\alpha$ 和 $\beta$。
  • 非线性谱减： 使用更复杂的减法规则（如基于后验信噪比的非线性函数）。
  • 多带谱减： 将频谱分成多个子带，在不同子带使用不同的减法参数。
  • 基于先验信噪比估计： 如最小均方误差（MMSE）谱减器（Ephraim & Malah, 1984），它利用先验信噪比来更平滑地控制谱减过程，显著减少音乐噪声。
  • 递归平均噪声估计： 在非语音段持续更新噪声估计，适应非平稳噪声。

优点

• 概念简单，计算复杂度相对较低。
• 对平稳或缓变噪声效果较好。

缺点

• 对非平稳噪声效果较差。
• 容易引入音乐噪声。
• 依赖于准确的语音活动检测（VAD）。
• 忽略了相位信息的影响。

基本流程

1. 信号转换：

   • 首先对时域信号进行短时傅里叶变换（STFT），将其转换到频域
   • 通常采用20-40ms的汉明窗（Hamming Window）进行分帧处理
   • 每帧之间有50%左右的重叠（overlap）

2. 噪声估计：

   • 在语音停顿段（即只有噪声的段落）计算噪声功率谱
   • 常用方法包括：前N帧平均法、最小值追踪法等
   • 噪声估计的准确性直接影响最终去噪效果

3. 谱减运算：

   • 从带噪语音功率谱中减去估计的噪声功率谱
   • 基本公式：|Y(ω)|² = |X(ω)|² - α|D(ω)|²
   • 其中α(0<α≤1)为过减因子，用于防止音乐噪声的产生

4. 后处理：

   • 对负值进行半波整流处理：|Y(ω)|² = max(|Y(ω)|², β|D(ω)|²)
   • β通常取0.01-0.1，称为谱下限因子
   • 最后通过逆STFT将信号转换回时域

应用场景

1. 语音通信系统（如VoIP、手机通话）
2. 语音识别前端处理
3. 助听器设备
4. 录音后期处理
5. 声学环境监测

实现示例（MATLAB伪代码）

% 1. 读取带噪语音
[noisy_speech, fs] = audioread('noisy.wav');

% 2. 分帧处理
frames = buffer(noisy_speech, frame_length, overlap);

% 3. 噪声估计（假设前5帧为纯噪声）
noise_spectrum = mean(abs(fft(frames(:,1:5))).^2, 2);

% 4. 谱减处理
for i = 1:size(frames,2)
    speech_fft = fft(frames(:,i));
    enhanced_mag = sqrt(max(abs(speech_fft).^2 - alpha*noise_spectrum, beta*noise_spectrum));
    enhanced_fft = enhanced_mag .* exp(1i*angle(speech_fft));
    enhanced_frames(:,i) = real(ifft(enhanced_fft));
end

% 5. 合成输出信号
enhanced_speech = overlap_add(enhanced_frames, overlap);

数据集地址

import os
import librosa
import librosa.display
import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf
from glob import glob

import warnings
from matplotlib import rcParams

rcParams['font.family'] = 'simhei'  # 解决图像的中文（黑体）显示问题
rcParams['axes.unicode_minus'] = False  # 解决负数坐标显示问题
warnings.filterwarnings('ignore')



# ---------------------- 1. NOIZEUS数据集加载工具 ----------------------
def load_noizeus_audio(noizeus_root, clean_speech_id="sp01", noise_type="babble", snr=0):
    """
    加载NOIZEUS数据集的纯净语音和对应带噪语音
    :param noizeus_root: NOIZEUS数据集根目录（包含clean/和noisy/文件夹）
    :param clean_speech_id: 纯净语音ID（如sp01~sp30）
    :param noise_type: 噪声类型（babble/car/white/office/restaurant/street/train）
    :param snr: 信噪比（0/5/10/15）
    :return: 纯净语音y_clean、带噪语音y_noisy、采样率sr
    """
    # 加载纯净语音
    clean_path = os.path.join(noizeus_root, f"{clean_speech_id}.wav")
    y_clean, sr = librosa.load(clean_path, sr=16000, mono=True)
    
    # 加载带噪语音（NOIZEUS文件名规则：spXX_噪声类型_snrX.wav）
    noisy_filename = f"{clean_speech_id}_{noise_type}_sn{snr}.wav"
    noisy_path = os.path.join(noizeus_root, noisy_filename)
    y_noisy, _ = librosa.load(noisy_path, sr=16000, mono=True)
    
    # 归一化
    y_clean = y_clean / np.max(np.abs(y_clean))
    y_noisy = y_noisy / np.max(np.abs(y_noisy))
    
    return y_clean, y_noisy, sr

def extract_noise_from_noizeus(y_noisy, sr, noise_duration=0.5):
    """
    从NOIZEUS带噪语音提取噪声样本（NOIZEUS开头约0.5秒为纯噪声）
    :param y_noisy: NOIZEUS带噪语音
    :param sr: 采样率（固定16kHz）
    :param noise_duration: 提取噪声时长（NOIZEUS推荐0.5秒）
    :return: 噪声样本y_noise
    """
    noise_samples = int(sr * noise_duration)
    y_noise = y_noisy[:noise_samples]
    return y_noise

# ---------------------- 2. NOIZEUS专属谱减法（16kHz最优参数） ----------------------
def spectral_subtraction_noizeus(y_noisy, y_noise, sr=16000):
    """
    适配NOIZEUS数据集的谱减法（16kHz最优参数）
    :param y_noisy: NOIZEUS带噪语音
    :param y_noise: 从带噪语音提取的噪声样本
    :param sr: 采样率（固定16kHz）
    :return: 降噪语音y_enhanced
    """
    # NOIZEUS最优参数（16kHz）
    n_fft = 512        # 32ms帧长，频率分辨率31.25Hz
    hop_length = 128   # 75%重叠率
    alpha = 1.05        # NOIZEUS噪声适配的过减系数
    
    # 1. 噪声STFT与平均幅度谱估计
    stft_noise = librosa.stft(y_noise, n_fft=n_fft, hop_length=hop_length)
    mag_noise_avg = np.mean(np.abs(stft_noise), axis=1, keepdims=True)
    
    # 2. 带噪语音STFT
    stft_noisy = librosa.stft(y_noisy, n_fft=n_fft, hop_length=hop_length)
    mag_noisy = np.abs(stft_noisy)
    phase_noisy = np.angle(stft_noisy)
    
    # 3. 谱减（加入过减因子和最小值补偿）
    beta = 1e-6
    mag_enhanced = np.maximum(mag_noisy - alpha * mag_noise_avg, beta)
    
    # 4. 逆STFT恢复语音
    stft_enhanced = mag_enhanced * np.exp(1j * phase_noisy)
    y_enhanced = librosa.istft(stft_enhanced, hop_length=hop_length)
    
    # 对齐长度+归一化
    y_enhanced = y_enhanced[:len(y_noisy)]
    y_enhanced = y_enhanced / np.max(np.abs(y_enhanced))
    
    return y_enhanced

# ---------------------- 3. 效果评估（SNR+STOI） ----------------------
def calculate_snr(clean, enhanced):
    """计算信噪比（SNR），评估降噪效果"""
    clean = clean[:len(enhanced)]
    noise = enhanced - clean
    snr = 10 * np.log10(np.sum(clean**2) / (np.sum(noise**2) + 1e-8))
    return snr

def calculate_stoi(clean, enhanced, sr=16000):
    """计算短时客观可懂度（STOI），需安装pystoi：pip install pystoi"""
    try:
        from pystoi import stoi
        clean = clean[:len(enhanced)]
        return stoi(clean, enhanced, sr, extended=False)
    except ImportError:
        print("未安装pystoi，跳过STOI计算（执行：pip install pystoi）")
        return None

# ---------------------- 4. 主流程：NOIZEUS谱减法实验 ----------------------
if __name__ == "__main__":
    # ========== 配置NOIZEUS参数 ==========
    NOIZEUS_ROOT = "./ori"  # 替换为你的NOIZEUS数据集根目录
    CLEAN_ID = "sp01"                  # 测试的纯净语音ID
    NOISE_TYPE = "car"              # 测试的噪声类型
    SNR_LEVEL = 0                      # 测试的信噪比（0dB难度最大）
    
    # ========== 1. 加载NOIZEUS数据 ==========
    y_clean, y_noisy, sr = load_noizeus_audio(
        noizeus_root=NOIZEUS_ROOT,
        clean_speech_id=CLEAN_ID,
        noise_type=NOISE_TYPE,
        snr=SNR_LEVEL
    )
    print(f"加载完成：纯净语音长度={len(y_clean)/sr:.2f}秒，采样率={sr}Hz")
    
    # ========== 2. 提取噪声样本 ==========
    y_noise = extract_noise_from_noizeus(y_noisy, sr, noise_duration=0.5)
    
    # ========== 3. 执行谱减法 ==========
    y_enhanced = spectral_subtraction_noizeus(y_noisy, y_noise, sr)
    
    # ========== 4. 保存结果 ==========
    save_dir = "noizeus_results"
    os.makedirs(save_dir, exist_ok=True)
    sf.write(os.path.join(save_dir, f"{CLEAN_ID}_{NOISE_TYPE}_snr{SNR_LEVEL}_clean.wav"), y_clean, sr)
    sf.write(os.path.join(save_dir, f"{CLEAN_ID}_{NOISE_TYPE}_snr{SNR_LEVEL}_noisy.wav"), y_noisy, sr)
    sf.write(os.path.join(save_dir, f"{CLEAN_ID}_{NOISE_TYPE}_snr{SNR_LEVEL}_enhanced.wav"), y_enhanced, sr)
    print(f"结果已保存到：{save_dir}")
    
    # ========== 5. 效果评估 ==========
    snr_noisy = calculate_snr(y_clean, y_noisy)
    snr_enhanced = calculate_snr(y_clean, y_enhanced)
    stoi_noisy = calculate_stoi(y_clean, y_noisy, sr)
    stoi_enhanced = calculate_stoi(y_clean, y_enhanced, sr)
    
    print("\n===== 降噪效果评估 =====")
    print(f"带噪语音SNR：{snr_noisy:.2f} dB")
    print(f"降噪后SNR：{snr_enhanced:.2f} dB")
    if stoi_noisy:
        print(f"带噪语音STOI：{stoi_noisy:.4f}")
        print(f"降噪后STOI：{stoi_enhanced:.4f}")
    
    # ========== 6. 可视化对比 ==========
    plt.figure(figsize=(15, 10))
    
    # 波形对比
    plt.subplot(3, 1, 1)
    librosa.display.waveshow(y_clean, sr=sr, color="#2E86AB")
    plt.title(f"纯净语音（{CLEAN_ID}）", fontsize=12)
    plt.ylabel("幅值")
    
    plt.subplot(3, 1, 2)
    librosa.display.waveshow(y_noisy, sr=sr, color="#A23B72")
    plt.title(f"带噪语音（{NOISE_TYPE}噪声，SNR={SNR_LEVEL}dB）", fontsize=12)
    plt.ylabel("幅值")
    
    plt.subplot(3, 1, 3)
    librosa.display.waveshow(y_enhanced, sr=sr, color="#F18F01")
    plt.title(f"谱减法降噪后语音", fontsize=12)
    plt.xlabel("时间 (s)")
    plt.ylabel("幅值")
    
    plt.tight_layout()
    plt.savefig(os.path.join(save_dir, f"{CLEAN_ID}_{NOISE_TYPE}_snr{SNR_LEVEL}_waveform.png"))
    plt.show()
    
    # 频谱对比（可选）
    plt.figure(figsize=(15, 5))
    n_fft = 512
    hop_length = 128
    
    # 带噪语音频谱
    plt.subplot(1, 2, 1)
    stft_noisy = librosa.stft(y_noisy, n_fft=n_fft, hop_length=hop_length)
    mag_noisy_db = librosa.amplitude_to_db(np.abs(stft_noisy), ref=np.max)
    librosa.display.specshow(mag_noisy_db, sr=sr, hop_length=hop_length, x_axis='time', y_axis='hz')
    plt.colorbar(format='%+2.0f dB')
    plt.title(f"带噪语音频谱（{NOISE_TYPE}，SNR={SNR_LEVEL}dB）", fontsize=12)
    
    # 降噪后频谱
    plt.subplot(1, 2, 2)
    stft_enhanced = librosa.stft(y_enhanced, n_fft=n_fft, hop_length=hop_length)
    mag_enhanced_db = librosa.amplitude_to_db(np.abs(stft_enhanced), ref=np.max)
    librosa.display.specshow(mag_enhanced_db, sr=sr, hop_length=hop_length, x_axis='time', y_axis='hz')
    plt.colorbar(format='%+2.0f dB')
    plt.title("谱减法降噪后频谱", fontsize=12)
    
    plt.tight_layout()
    plt.savefig(os.path.join(save_dir, f"{CLEAN_ID}_{NOISE_TYPE}_snr{SNR_LEVEL}_spectrogram.png"))
    plt.show()

相关论文数据集与实验结果

谱减法及其改进算法被广泛应用于各种语音数据集上进行评估。以下是一些经典论文中提到的数据集和典型实验结果：

1. 论文示例 1： Boll, S. F. (1979). Suppression of Acoustic Noise in Speech Using Spectral Subtraction. IEEE Transactions on Acoustics, Speech, and Signal Processing, 27(2), 113–120.

   • 数据集： 通常使用实验室录制的纯净语音叠加人工添加的噪声（如白噪声、粉红噪声、办公室噪声）或真实环境噪声录音。具体数据集名称文中可能未明确给出（当时标准数据集较少），但方法本身是通用的。
   • 实验结果： 论文展示了谱减法在信噪比（SNR）提升方面的有效性。例如，在输入SNR为5dB的白噪声环境下，谱减法处理后输出SNR可能提升到10dB左右（具体数值取决于参数设置和噪声类型）。主观听感报告噪声被有效抑制，但存在可感知的失真和残留噪声（音乐噪声的雏形）。

2. 论文示例 2： Ephraim, Y., & Malah, D. (1984). Speech Enhancement Using a Minimum Mean-Square Error Short-Time Spectral Amplitude Estimator. IEEE Transactions on Acoustics, Speech, and Signal Processing, 32(6), 1109–1121.

   • 数据集： 该篇提出MMSE-STSA算法的经典论文也使用了类似的人工加噪方式。后续研究广泛使用标准数据集进行评测，例如：
     • TIMIT: 包含多种方言的纯净语音数据库，常用于叠加噪声。
     • NOISEX-92: 包含多种类型的环境噪声（如白噪、粉红噪、babble噪、工厂噪、汽车噪）。
   • 实验结果： MMSE-STSA（常被视为一种更先进的谱减框架）相对于基本谱减法有显著改进。论文报告了在多种噪声和SNR条件下，MMSE-STSA在客观指标（如分段SNR, SegSNR）上的提升远优于基本谱减法。例如，在5dB白噪声下，基本谱减法SegSNR提升可能只有2-3dB，而MMSE-STSA可能达到5-6dB的提升。更重要的是，主观听感测试（如MOS - Mean Opinion Score）显示MMSE-STSA产生的音乐噪声远少于基本谱减法，语音自然度和清晰度更高。例如，MOS得分（1-5分，越高越好）可能从基本谱减法的2.5左右提高到MMSE-STSA的3.5以上（具体取决于噪声条件和SNR）。

总结

谱减法是语音增强领域的基石之一。虽然基本形式存在音乐噪声等问题，但其思想启发了大量后续研究（如MMSE系列算法）。在平稳噪声环境下，即使是基本谱减法也能有效提升SNR。改进的谱减法（如MMSE-STSA）在客观指标和主观听感上都有显著提升，并被用作后续深度学习语音增强方法的基准对比算法之一。评估这些算法常用的数据集包括TIMIT（纯净语音）和NOISEX-92（噪声），实验结果通常报告客观SNR提升和主观MOS评分。