本文简要探讨了语音信号的基本特性，并介绍了谱减法的基本原理，同时基于Matlab实现了原始谱减法及其一种改进版。最后给出音频波形、频谱，从主观的人耳和客观的评判标准来量化所达到的效果。

1 语音信号

1.1 宏观的非平稳性与微观的平稳性

在现实世界中，我们获得的所有信号波形都是其对应分布的一组观测值。平稳信号指的是其分布及分布的参数不会发生变化，非平稳信号则是其分布或分布的参数发生变化。对于语音信号来说，由于信号直接由说话者控制，其分布和分布的参数都是不确定的。以一个宏观的角度来看，语音信号具有非平稳性。而在 $10$ - $30ms$ 这个尺度，受人的发声器官所限定，语音信号是平稳的，即在微观下，语音信号具有平稳性。正是这种特性，导致了在语音信号处理过程中，短时分析成为一个重要的手段。

2 基本谱减法的原理

对噪声的基本假设

噪声是局部稳定的，即与语音信号一致，在微观尺度下具有平稳性。
噪声是加性噪声
噪声与纯净的语音信号是不相关的

基于以上的假设，我们如下定义

$x(t)$ 为纯净的语音信号
$n(t)$ 为噪声信号
$y(t)$ 为加噪后的语音信号

则有

$y(t) = x(t) + n(t)$

对加噪后的语音信号做傅里叶变换可得

$Y(\omega) = X(\omega) + N(\omega)$

可得加噪后的能量谱如下

$|Y(\omega)|^2 = |X(\omega)|^2 + |N(\omega)|^2 +X^*(\omega)N(\omega) + N^*(\omega)X(\omega)$

又因为噪声与纯净的语音信号不相关，则

$|Y(\omega)|^2 = |X(\omega)|^2 + |N(\omega)|^2$

即

$|X(\omega)|^2 = |Y(\omega)|^2 - |N(\omega)|^2$

亦由于语音信号的质量与其频率幅度谱密切相关，而与频率相位图并无太大关联。则，我们可以保留加噪后的语音信号的频率相位谱，并且估计出噪声的能量谱，就通过以上式子进行语音增强

3 论文复现

我们复现了M.Berouti等人的Enhancement of speech corrupted by acoustic noise中提出的一种改进版本的谱减法，并使用了Dong Wang等人发布的开源语言数据集THCHS-30

3.1 音频波形展示

时域信号幅度图

语谱图

3.2 代码

3.2.1 谱减法

function output = spectral_subtraction(input, wlen, inc, nlen, a, b)
%------------------------------------
%  input : input signal
%  wlen  : the size of each frame
%  inc   : the size of frame shift
%  nlen  : the size of noisy without speech
%  output: output signal
%------------------------------------
window_function = hamming(wlen);                       % set window function

y = enframe(input, window_function, inc)';             % subframe: y is a matrix with wlen×frame_n
frames_n = size(y, 2);                                 % the number of frames

y_fft = fft(y);                                        % fft
y_fft_amplitude = abs(y_fft);                          % amplitude
y_fft_energy = y_fft_amplitude.^2;                     % energy
y_fft_phase = angle(y_fft);                            % phase

noisy_energy = mean(y_fft_energy(:,1:nlen), 2);        % noisy energy is a vector with wlen×1

output_energy = ones(wlen, frames_n);
for ii = 1:frames_n
    for jj = 1:wlen
        if y_fft_energy(jj, ii) > a*noisy_energy(jj)
            output_energy(jj, ii) = y_fft_energy(jj, ii) - a*noisy_energy(jj);
        else
            output_energy(jj, ii) = b*y_fft_energy(jj, ii);
        end        
    end
end
output_amplitude = sqrt(output_energy);                % output amplitude

output = OverlapAdd2(output_amplitude, y_fft_phase, wlen, inc);

end

谱减法（Spectral Subtraction）