subspectrum.m

ＭATLAB代码,利用谱减法实现语音信号增强.经过测试,效果不错.

function enhancedsignal=subspectrum(noisysignal,samplefrequency)
x=noisysignal;
fs=samplefrequency;
nx=length(x);
enhanced_x=zeros(1,nx);
wingain=zeros(1,nx);

%分帧和加窗Framing and windowing
FrameLen=fix(0.025*fs);         %帧长the length of frame
overlap=FrameLen/2;
inc=FrameLen-overlap;           %帧移
x_frame=enframe(x,FrameLen,inc); %分帧framing
nf=size(x_frame,1);              % 帧数the number of frames
win=hamming(FrameLen)';
x_window=[];
for k=1:nf
    x_row=x_frame(k,:).*win;     % 加窗windowing
    x_window=[x_window;x_row];   
end

%对带噪语音进行DFT 
y=fft(x_window');
ymag = abs(y);			
yphase = angle(y);	

%初始噪声估计 Noise estimation
NNoise=20; %取噪音段（语音的初始段）帧数 get a segment of speech 
MN=mean(ymag(:,1:NNoise)')';
NoiseCounter=0;%连续噪声段长度 the length of noise 
SmoothFactor=9;%噪声平滑因子  smoothing factoor

%谱减实现 spectral subtraction 
for k=1:nf
%-----噪声估计更新 update of noise estimation
     if k<=NNoise 
        SpeechFlag=0;%非有声段 not speech
        NoiseCounter=NNoise;
    else 
           NoiseMargin=3;
           HangOver=8;
           SpectralDist= 20*(log10(ymag(:,k))-log10(MN));
           SpectralDist(find(SpectralDist<0))=0;
           Dist=mean(SpectralDist); 
           if (Dist < NoiseMargin) 
            NoiseFlag=1; 
            NoiseCounter=NoiseCounter+1;
           else
            NoiseFlag=0;
            NoiseCounter=0;
           end
           % 只检测大于一定长度的噪声段 % Detect noise only periods and attenuate the signal
           if (NoiseCounter > HangOver) 
           SpeechFlag=0;    
           else 
           SpeechFlag=1; 
           end 
     end
    
    if SpeechFlag==0 % 如果是噪声段% If not Speech Update Noise Parameters
        MN=(SmoothFactor*MN+ymag(:,k))/(SmoothFactor+1); %更新噪声均值
    end
%-----幅度谱减
 smag=ymag(:,k)-MN;
 smag(find(smag<0))=0;
%------重新生成谱
 spectrum= smag.*exp(j*yphase(:,k));
%------重新生成语音
 enhanced_x((inc*(k-1)+1):(inc*(k-1)+FrameLen))=enhanced_x((inc*(k-1)+1):(inc*(k-1)+FrameLen))+real(ifft(spectrum))';
%------窗增益
 wingain(1+inc*(k-1):FrameLen+inc*(k-1))=wingain(1+inc*(k-1):FrameLen+inc*(k-1))+win; 
end

% 除去Hamming窗引起的增益
for k=1:nx
    if  wingain(k)==0
        enhanced_x(k)=0;
    else
        enhanced_x(k)=enhanced_x(k)/wingain(k);
    end
end 
enhancedsignal=enhanced_x;