音频噪声抑制（7）：利用基本谱减法实现噪声抑制

高二的时候吧，我们拍过一部校园DV剧…《十六岁的伤悲》……

班上的大叔用Cool Edit Pro这个软件，给剧中的旁白去了噪声。

大致过程是这样的（Cool Edit Pro 去噪）

1. 在录音的时候，先空出几秒不说话，相当于是在录环境噪声；

2. 然后开始录配音演员的声音；

3. 用Cool选出前面那段噪声，分析一下噪声的参数；

4. 把整段音频选中，导入刚刚分析出来的参数实现降噪处理。

有点神奇。

7年过去了，我现在想，是不是大概靠的谱减法去噪？

当然，一般情况下，也没这么简单，毕竟降噪效果挺好的。

不过，千里之行始于足下，所以还是先学习基本的谱减法。

在一本书上，找到了详细的介绍。

以下的符号中，y是含噪信号，S是干净信号，N是加性噪声。









注意。

1. 在去噪过程中，对语音进行了加窗、分帧，因为有个窗长（window size）、步进长度（hop size），所以，在处理完之后，不能直接相加，而是应该采用重叠相加（overlapped add），重叠相加之前，还应该去掉因加窗造成的加窗增益（在以下代码中，还暂时没有考虑这一点）。一定一定要重叠相加！不然会出现致命的问题！

2. 在实际录音过程中，噪声也不会这么理想，就是加性噪声。万一是乘性噪声呢？是不是应该考虑取个对数，把乘法运算降级为加法运算，处理完以后再指数翻上去？

这种“降级”非常常见。如，利用log把乘变为加，利用fourier/laplace transform把卷积（积分）变为乘法，……

我用手机录音的时候，还是有一定的环境噪声的（笔记本散热，风扇转得很厉害，呼啦呼啦的）。

%% 谱减法去噪
% 内容：用最简单的谱减法，去除噪声
% 作者：qcy
% 版本：v1.0
% 【注意】谱减完以后，拼接在一起的时候，要去除重复的部分！
% 因为分帧时有可能是有重叠的！！！
% 时间：2016年11月3日19:40:35

clear; clc; close all;
%% 导入音频

filedir=[];                             % 指定文件路径
filename='静默检测.m4a';                % 指定文件名
fle=[filedir filename];                 % 构成路径和文件名的字符串
[xx,fs]=audioread(fle);                 % 读入数据文件
xx=xx-mean(xx);                         % 消除直流分量
x=xx/max(abs(xx));                      % 幅值归一化

%% 预先设置参数
IS = 1.6;      % 设置前导无话段长度 [s]
wlen = 200;    % 务必设置为偶数个点
inc = 80;
win = hamming(wlen); % 用矩形窗，也会有问题-_-!
% boxcar 矩形窗
% hamming 汉明窗

N=length(x);                            % 信号长度
time=(0:N-1)/fs;                        % 设置时间
overlap=wlen-inc;                       % 求重叠区长度
NIS=fix((IS*fs-wlen)/inc +1);           % 求前导无话段帧数
Nframe = floor( (length(x) - wlen) / inc) + 1; % 一共多少帧

%% 分帧
% 1) 测量环境噪声的平均能量谱（只需计算非负频率的部分）
k_pos_freq = wlen/2+1; % 非负频率范围
X_noise_engergy_sum = zeros(k_pos_freq,1);

for k = 1 : NIS
idx = (1:wlen) + (k-1) * inc;
x_temp = x(idx).*win;
X_temp = fft(x_temp);
X_noise_engergy_sum = X_noise_engergy_sum + abs(X_temp(1:k_pos_freq)).^2;
end
X_noise_engergy_avg = X_noise_engergy_sum / NIS;

% 2) 对之后语音每一帧进行谱减
% 3) 注意：重构的时候要重叠相加-_-!
a = 10;
b = 0.002;

Y_ = zeros(wlen,1);
x_subspec = [];

sig=zeros((Nframe-1)*inc+wlen,1);

for k = 1 : Nframe
idx = (1:wlen) + (k-1) * inc;
x_temp = x(idx).*win;
X_temp = fft(x_temp);
phase_k = angle(X_temp(1:k_pos_freq)); % 1. 先保留信号的相位
X_energy = abs(X_temp(1:k_pos_freq)).^2;

X_subspec_energy_pos_freq = zeros(k_pos_freq,1); % 谱减以后，非负频率分量的能量
for m = 1:k_pos_freq % 对每一个频点进行谱减
if  X_energy(m) >= a * X_noise_engergy_avg(m)
X_subspec_energy_pos_freq(m) = X_energy(m) - a * X_noise_engergy_avg(m);
else
X_subspec_energy_pos_freq(m) = b * X_noise_engergy_avg(m);
end
end

X_subspec_abs_pos_freq = sqrt(X_subspec_energy_pos_freq);
X_subspec_pos_freq = X_subspec_abs_pos_freq .* exp(1j*phase_k); % A = |A|exp(1j*angle(A))
% 构造Hermitian对称的序列
X_subspec_at_frame_k = [X_subspec_pos_freq ; conj(X_subspec_pos_freq(end-1:-1:2))];
x_subspec_at_frame_k = real(ifft(X_subspec_at_frame_k));

%     x_subspec = [x_subspec ; x_subspec_at_frame_k ]; % 直接接上是不行的

% 重叠相加 -_-! 为什么要重叠相加!!才不会有嘶哑的声音!!!!
% 由频域重构 -_-! 难道是说，加窗以后变到频域，再变回来，重构，需要重叠相加抵消xx效应?

start=(k-1)*inc+1;
sig(start:start+wlen-1)=sig(start:start+wlen-1) + x_subspec_at_frame_k; % 重叠相加
end

%%
figure;
subplot(211);
plot(x);
title('去噪前');
subplot(212);
plot(sig);
title('去噪后');
sound(sig,fs);