specsubtd_boll79.m

基于谱减法的语音增强在MATLAB实现源程序

% Spectral Subtraction in time domain ================
%
% This is an implementation of the Spectral Subtraction method described in 
% Boll 79.
%
% Usage:
%   x_hat = SpecSubTD_Boll79(y);
% Example:
%   y = wavread_aurora(/aurora/speechdata/TESTA/N1_SNR5/FAK_1B.08);
%   x_hat = SpecSubTD_Boll79(y);
%
% input:
%   y: a Nx1 matrix of samples (assume 8kHz)
%
% output:
%   x_hat: a Nx1 matrix of samples
%
% Updates and other algorithms available at: http://www.RobustASR.net
%
% Released under GNU GPL license
% Copyright 2001, Trausti Kristjansson
% ===================================================

% Overview of algorithm:
%
% take fft of overlapping frames of time domain signal (assume 8kHz)
%   1) Smooth data with simple 3 tap filter
%      Estimate avarage noise power
%      Find residual noise power
%   2) Spectral subtraction
%   3) Half wave rectification
%   4) Residual Noise Reduction
%   5) Additional Signal Attenuation During Nonspeech Activity
% take ifft and do overlap add to convert back into time domain

function [X_hat]=SpecSubTD_Boll79(y);
[y]=prosnr(y,0);
n_fl=ones(256,1);
% Take the FFT of overlapping frames
wind=hanning(256);
i=1;
for t=1:128:(length(y)-256);
    sp(:,i)=fft(wind.*y(t:(t+255)));
    abs_sp(:,i)=abs(sp(:,i));
    angle_sp(:,i)=angle(sp(:,i));
    i=i+1;
end;

Y_FB=abs_sp; 

% 1) Magnitude Averaging =========================== 
% simple 3 tap lowpass filtering
Y_FBav=Y_FB;
for i=2:(size(Y_FB,2)-1)
    Y_FBav(:,i)=mean(Y_FB(:,(i-1):(i+1)),2);
end;
Y_FB=Y_FBav;

% use first 20 frames to estimate the noise mean

N_FB_mean=mean(Y_FB(:,1:20),2);

% Find maximum deviation for use in Residual Noise Reduction step

N_R=max(Y_FB(:,1:20)-N_FB_mean*ones(1,20),[],2);

% 2) Spectral Subtraction ==========================
X_hat_FB=Y_FB-(N_FB_mean-n_fl)*ones(1,size(Y_FB,2));

% 3) Half wave rectification =======================
% set values under a threshold to threshold
n_flooro=n_fl*ones(1,size(Y_FB,2));
[I,J]=find(X_hat_FB<n_flooro);
X_hat_FB(sub2ind(size(X_hat_FB),I,J))=n_flooro(sub2ind(size(X_hat_FB),I,J));


% 4) Residual Noise Reduction ====================== 
% A type of median filtering:
% if noise is less than N_R, then use min value over adjacent frames
X_hat_FBrn=X_hat_FB;
for i=2:(size(Y_FB,2)-1)
    I=find( X_hat_FB(:,i)<(N_R + n_fl) );
    X_hat_FBrn(I,i)=min(X_hat_FB(I,(i-1:i+1)),[],2);
end;
X_hat_FB=X_hat_FBrn;

% 5) Additional Signal Attenuation During Nonspeech Activity
% This step classifies frames into speech/non-speech,
% and replaces non-speech frames by attenuaated versions of the
% original frame. The replaced fames have power equivalent to
% the noise floor.
% find attenuation constant
c=10^(-30/20);

% find non-speech frames
T=20*log10(mean(X_hat_FB./(N_FB_mean*ones(1,size(X_hat_FB,2)))));
noiseFrames = find(T<30);

% replace them with attenuated noisy frames.
X_hat_FB(:,noiseFrames)=c*Y_FB(:,noiseFrames);

% convert back to time domain, by overlap-add

X_hat=zeros(1,length(y));
j=sqrt(-1);
i=1;
for t=1:128:(length(y)-256);
    X_hat(:,(t:t+255))=X_hat(:,(t:t+255))+real(ifft(X_hat_FB(:,i).*exp(j*angle_sp(:,i))))';
    i=i+1;
end;
%xs=x_hat(1,:);
%for i=2:4%((length(y)-256)/128)
%    xs=[xs x_hat(i,129:256)];%c=[c a(3,2:3)]
%end;
subplot(211)
plot(y)
subplot(212)
plot(X_hat)