subspec_magic.m

基于谱减法的语音增强程序

function subspec_magic(a,b);
%经典谱减法  语音增强程序
% 07.12.20  magic

[filename,pathname]=uigetfile('*.wav','请选择语音文件：');
[wavin,fs,nbits]=wavread([pathname filename]);

%基音周期最大为20ms，为使ifft还原后语音失真尽量小，帧长至少要为基音周期的2倍
%根据fs选择帧长：
switch fs
    case 8000
        frame_len=320;step_len=160;
    case 10000
        frame_len=400;step_len=200;
    case 12000
        frame_len=480;step_len=240;
    case 16000
        frame_len=640;step_len=320;
    case 44100
        frame_len=1800;step_len=900;
    otherwise
        frame_len=1800;step_len=900;
end;
inframe=(ENFRAME(wavin,frame_len,step_len))';   %分帧
frame_num=size(inframe,2);          %求帧数
window=hamming(frame_len);          %定义汉明窗
% window=boxcar(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;

%相邻3帧平滑
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;

%取前20帧，求平均，做为噪声幅值
abs_noise=mean(abs_frame(:,1:20),2);

%对前20帧求最大残余噪声
res_noise=max(abs_frame(:,1:20)-abs_noise*ones(1,20),[],2);

%谱减运算，可适当调节a,b的值来获取合适的结果
%(sub_frame=abs_frame^a-b*abs_noise^a)
if a==1
     for i=1:frame_num;
        sub_frame(:,i)=abs_frame(:,i)-b*abs_noise;
     end;
else for i=1:frame_num;
        sub_frame_p(:,i)=abs_frame(:,i).^a-b*(abs_noise.^a);
        sub_frame(:,i)=sub_frame_p(:,i).^(1/a);
    end;
end;

zero_frame=zeros(frame_len,frame_num);       %定义与sub_frame同维数的零距阵
sub_frame_z=sub_frame;

%谱减结果修正：将sub_frame距阵中小于0的点置0
[I,J]=find(sub_frame_z<zero_frame);
sub_frame_z(sub2ind(size(sub_frame_z),I,J))=zero_frame(sub2ind(size(sub_frame_z),I,J));
sub_frame=sub_frame_z;

%残余噪声消去：将sub_frame中小于res_noise的点用相邻3中的最小值来代替
for i=2:(frame_num-1);
    I=find(sub_frame(:,i)<res_noise);
    sub_frame(I,i)=min(sub_frame_z(I,(i-1):(i+1)),[],2);
end;

% 非语音帧衰减
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);
% noise_frame=find(T<-20);
sub_frame(:,noise_frame)=c*sub_frame(:,noise_frame);

%将处理结果还原至时域
wavout=zeros(1,(frame_num-1)*step_len+frame_len);
j=sqrt(-1);
% for i=1:frame_len;
%     wavout(((i-1)*step_len+1):((i-1)*step_len+frame_len))=...
%     wavout(((i-1)*step_len+1):((i-1)*step_len+frame_len))+real(ifft(sub_frame(:,i).*(exp^(j*ang_frame(:,i))))';
% end;
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'文件
wavwrite(wavout,fs,nbits,['2ss_' num2str(frame_len) 'f_' num2str(a) num2str(b) '_' filename]);

%将处理前后的结果进行作图比较
subplot(2,1,1);
plot(wavin);grid on;
subplot(2,1,2);
plot(wavout);grid on;