4. 谱减法混响抑制
任务说明
本案例实现基于谱减法的混响抑制算法。输入为带混响的语音信号（由干净语音与人工混响冲激响应卷积生成），通过短时傅里叶变换（STFT）将信号变换到频域，估计噪声/混响功率谱，执行谱减法得到增强频谱，最后通过逆STFT恢复时域信号。重点比较不同过减因子（over-subtraction factor）α对混响抑制效果的影响，以频谱图与信混比（SRR）提升为判定依据。
算法原理与实现
谱减法的核心思想是在频域中将估计的混响功率谱从观测信号功率谱中减去。设带混响信号为 y(n)，其STFT为 Y(k,l)，其中 k 为频率索引，l 为帧索引。增强信号频谱估计为：
|X̂(k,l)|² = |Y(k,l)|² − α·|D(k,l)|²
其中 |D(k,l)|² 为混响功率谱估计，α 为过减因子。当 α 增大时，混响被更激进地去除，但可能引入音乐噪声或语音失真。本实现中，混响功率谱采用带混响信号前若干帧（假设混响尾音主要在起始段）的统计平均作为估计，并设置谱下限防止负值。
import numpy as np
import soundfile as sf
from scipy import signal
import matplotlib.pyplot as plt
def stft(x, win_len=512, hop_len=128):
"""短时傅里叶变换，返回复数频谱矩阵 (帧数, 频点数)"""
n_frames = 1 + (len(x) - win_len) // hop_len
win = np.hanning(win_len)
X = np.zeros((n_frames, win_len//2 + 1), dtype=complex)
for i in range(n_frames):
start = i * hop_len
seg = x[start:start+win_len] * win
X[i, :] = np.fft.rfft(seg)
return X
def istft(X, win_len=512, hop_len=128):
"""逆短时傅里叶变换，重叠相加恢复时域信号"""
n_frames, freq_bins = X.shape
x = np.zeros((n_frames - 1) * hop_len + win_len)
win = np.hanning(win_len)
for i in range(n_frames):
seg = np.fft.irfft(X[i, :], n=win_len)
start = i * hop_len
x[start:start+win_len] += seg * win
# 归一化（窗函数重叠相加的补偿）
norm = np.zeros_like(x)
for i in range(n_frames):
start = i * hop_len
norm[start:start+win_len] += win**2
norm[norm < 1e-8] = 1.0
return x / norm
def spectral_subtraction(y, alpha, n_fft=512, hop=128, n_est_frames=10):
"""谱减法混响抑制
y: 带混响信号
alpha: 过减因子
n_est_frames: 用于估计混响功率谱的起始帧数
"""
Y = stft(y, n_fft, hop)
mag_Y = np.abs(Y)
phase_Y = np.angle(Y)
# 用前n_est_frames帧的平均功率作为混响功率谱估计
noise_pow = np.mean(mag_Y[:n_est_frames, :]**2, axis=0)
# 谱减
mag_X2 = mag_Y**2 - alpha * noise_pow
mag_X2 = np.maximum(mag_X2, 0.01 * noise_pow) # 谱下限防止负值
mag_X = np.sqrt(mag_X2)
X_enh = mag_X * np.exp(1j * phase_Y)
return istft(X_enh, n_fft, hop)
def compute_srr(clean, enhanced):
"""计算信混比提升（简化：用能量比近似）"""
# 将增强信号与干净信号对齐后计算误差
min_len = min(len(clean), len(enhanced))
c = clean[:min_len]
e = enhanced[:min_len]
# 简单对齐：找最大相关延迟
corr = np.correlate(e, c, mode='full')
delay = np.argmax(corr) - (len(c) - 1)
if delay > 0:
e_aligned = e[delay:]
c_aligned = c[:len(e_aligned)]
else:
e_aligned = e[:len(c)+delay]
c_aligned = c[-delay:]
min_l = min(len(c_aligned), len(e_aligned))
c_aligned = c_aligned[:min_l]
e_aligned = e_aligned[:min_l]
err = c_aligned - e_aligned
srr = 10 * np.log10(np.sum(c_aligned**2) / (np.sum(err**2) + 1e-10))
return srr
# ========== 主流程 ==========
# 1. 生成测试数据：干净语音 + 人工混响
fs = 16000
t = np.arange(0, 2.0, 1/fs)
clean_speech = 0.5 * np.sin(2*np.pi*220*t) * np.exp(-t*2) # 衰减正弦模拟语音
clean_speech += 0.3 * np.sin(2*np.pi*440*t) * np.exp(-t*3)
# 人工混响冲激响应（指数衰减）
rir_len = 4000
rir = np.exp(-np.arange(rir_len)/2000) * np.random.randn(rir_len)
rir /= np.sqrt(np.sum(rir**2))
reverb_speech = signal.fftconvolve(clean_speech, rir)[:len(clean_speech)]
# 混合（混响能量比例约50%）
y = clean_speech + 0.8 * reverb_speech
# 2. 不同过减因子比较
alphas = [0.5, 1.0, 2.0, 4.0]
results = {}
for a in alphas:
enhanced = spectral_subtraction(y, alpha=a)
results[a] = enhanced
srr = compute_srr(clean_speech, enhanced)
print(f"alpha={a:.1f}, SRR提升={srr:.2f} dB")
# 3. 可视化对比（保存图片）
plt.figure(figsize=(12, 8))
plt.subplot(3, 2, 1)
plt.specgram(clean_speech, Fs=fs, NFFT=512, noverlap=384)
plt.title("Clean Speech")
plt.subplot(3, 2, 2)
plt.specgram(y, Fs=fs, NFFT=512, noverlap=384)
plt.title("Reverberant Speech")
for idx, a in enumerate(alphas):
plt.subplot(3, 2, idx+3)
plt.specgram(results[a], Fs=fs, NFFT=512, noverlap=384)
plt.title(f"Enhanced (alpha={a})")
plt.tight_layout()
plt.savefig("spectral_subtraction_result.png", dpi=100)
print("频谱图已保存至 spectral_subtraction_result.png")运行与验证
依赖安装：
pip install numpy scipy soundfile matplotlib运行方式：将上述代码保存为 spec_subtract.py，在Python环境中执行 python spec_subtract.py。
验证方法：
- 观察控制台输出的SRR提升值：随着α从0.5增大到4.0，SRR应呈现先升后降的趋势（α过小抑制不足，α过大引入失真）。
- 查看生成的
spectral_subtraction_result.png：带混响语音的频谱图呈现明显的“拖尾”模糊现象；增强后频谱中拖尾被削弱，但α较大时可能出现零星噪声点（音乐噪声）。 - 判定标准：α=1.0~2.0区间内SRR提升达到峰值（通常比α=0.5高3~6dB），且频谱图清晰度明显改善，即视为算法有效。