untitled.m

提供了集中语音去噪的matlab代码

%本程序应用多窗谱法估计的语音信号功率谱密度（PSD）来进行谱减语音增强

clear;
a=2;      %过减因子
b=0.01;     %增益补偿因子
c=0;        %c=0时，不对增益矩阵进行开方，c=1时，进行开方运算
n=0.04;
%读取语音文件---------------------------------------------------------------
[filename,pathname]=uigetfile('*.wav','请选择语音文件：');
[wavin,fs,nbits]=wavread([pathname filename]);
wav_length=length(wavin);
[speech,fs,nbits]=wavread('speech_clean.wav');         % 读入数据 

%基音周期最大为20ms，为使ifft还原后语音失真尽量小，帧长至少要为基音周期的2倍
%根据fs选择帧长：
frame_len=256;step_len=128;
inframe=(enframe(wavin,frame_len,step_len))';   %分帧

frame_num=size(inframe,2);          %求帧数
window=hamming(frame_len);          %定义汉明窗

%分别对每帧fft，求幅值，求相角-----------------------------------------------
for i=1:frame_num;
    fft_frame(:,i)=fft(window.*inframe(:,i));
    abs_frame(:,i)=abs(fft_frame(:,i));
    ang_frame(:,i)=angle(fft_frame(:,i));
end;

%每相邻三帧平滑-------------------------------------------------------------
abs_frame_f=abs_frame;
for i=2:frame_num-1;
    abs_frame_f(:,i)=mean(abs_frame(:,(i-1):(i+1)),2);
end;
abs_frame=abs_frame_f;

%求增益矩阵-----------------------------------------------------------------
%矩阵中每一元素为：
%g(k)=(Py(k)-a*Pn(k))/Py(k)
%Py和Pn分别为带噪语音和噪声的功率谱估计，都用MATLAB中自带的pmtm函数来估计
%可根据需要调节a的大小，来得到更好的效果

%用多窗谱法法对每一帧数据进行功率谱估计
for i=1:frame_num;
    per_PSD(:,i)=pmtm(inframe(:,i),3,frame_len,'twosided');
end;

%对功率谱的每相邻三帧进行平滑
per_PSD_f=per_PSD;
for i=2:frame_num-1;
    per_PSD_f(:,i)=mean(per_PSD(:,(i-1):(i+1)),2);
end;
per_PSD=per_PSD_f;

%取前20帧作为噪声帧，取其平均作为噪声的功率谱估计
% noisy=inframe(1:20);
noise_PSD=mean(per_PSD(:,1:20),2);
abs_noise=mean(per_PSD(:,1:20),2).^0.5;
%求增益矩阵
for k=1:frame_num;
    g(:,k)=(per_PSD(:,k)-a*noise_PSD)./per_PSD(:,k);
end;

%求增益补偿阈值，凡是小于该阈值的增益系数军用阈值来代替，这样可减少音乐噪声
spec_floor=b*noise_PSD./per_PSD(:,k);
spec_floor=spec_floor(:,ones(1,frame_num));
[I,J]=find(g<spec_floor);
gf=g;
gf(sub2ind(size(gf),I,J))=spec_floor(sub2ind(size(gf),I,J));
if c==0;
    g=gf;
else g=gf.^0.5;
end;
% g=gf.^0.5;
%谱减----------------------------------------------------------------------
sub_frame=g.*abs_frame;

% 非语音帧衰减------------------------------------------------------------
T=20*log10(mean(sub_frame./(abs_noise*ones(1,frame_num))));
T_noise=mean(T(:,1:20),2);
c=10^(-2/3);            %衰减系数为 10^(-1.5)
noise_frame=find(T<T_noise);
sub_frame(:,noise_frame)=c*sub_frame(:,noise_frame);

%将语音信号还原至时域-------------------------------------------------------
wavout=zeros(1,(frame_num-1)*step_len+frame_len);
j=sqrt(-1);
i=1;
for t=1:step_len:((frame_num-1)*step_len+1);
    wavout(:,t:(t+frame_len-1))=wavout(:,t:(t+frame_len-1))+real(ifft(sub_frame(:,i).*exp(j*ang_frame(:,i))))';
    i=i+1;
end;

%将处理结果输出为'wav'文件--------------------------------------------------
if c==0;
    wavwrite(wavout,fs,nbits,[num2str(frame_len) 'fnm_' num2str(a) '_' num2str(b) '_' filename]);
else wavwrite(wavout,fs,nbits,['m_' num2str(a) '_' num2str(b) '_' filename]);
end;

SNR1 = 10*log10(var(speech)/var(wavin-speech));               %加噪语音信噪比

waveout=[wavout  zeros(1,(length(wavin)-length(wavout')))];

SNR2 = 10*log10(var(speech')/var(waveout'-speech));     %增强语音信噪比

% %将处理前后的结果进行作图比较
% figure(1)
% subplot(3,1,1);
% plot(wavin);grid on;
% title(['Noise Added (SNR=',num2str(SNR1),'dB)']);
% axis([1 wav_length -1 1]);
% subplot(3,1,2);
% plot(wavout);grid on;
% title(['Improve Voice (SNR=',num2str(SNR2),'dB)']);
% axis([1 wav_length -1 1]);
% subplot(313);
% plot(speech);grid on;
% title(['Original Voice (n=',num2str(n),')' ]);
% axis([1 wav_length -1 1]);
% %将处理前后的结果进行作图比较
% figure(2)
% specgram(speech');title(['Original Voice (n=',num2str(n),')' ]);
% figure(3)
% specgram(wavin);  title(['Noise Added (SNR=',num2str(SNR1),'dB)']);
% figure(4)
% specgram(wavout); title(['Improved Voice (SNR=',num2str(SNR2),'dB ¤SNR=',num2str(SNR2-SNR1),'dB)']);
% 
% pause;
% soundsc(wavin,fs,nbits);
% pause;
% soundsc(wavout,fs,nbits);
% pause;
% soundsc(speech,fs,nbits);
%