evrc_5.m

matlab仿真通过的降噪程序

%%%%%%%  add:  update  see EVRC IS 127  %%%%%%%

clear all;
[y1,fs,bits]=wavread('E:\noise enhancing\wav\5.wav');
y1=y1/max(abs(y1));%归一化
%wavwrite(y1,8000,8,'f:\wav\3music.wav');
figure(1);
plot(y1);

[noise,fs1,bits1]=wavread('E:\noise enhancing\wav\5_noise.wav');
y2=mixsig(y1,noise,5);
y2=y2/max(abs(y2));%归一化
wavwrite(y2,8000,8,'E:\noise enhancing\wav\mymasking_s&w(10).wav');%0db带噪信号
figure(2);
plot(y2);

% window=chebwin(35,80);
% b=fir1(34,0.015,'high',window);
% y=filter(b,1,y2);
[b,a]=butter(6,120/4000,'high');
y=filter(b,a,y2);
figure(3);
plot(y);
wavwrite(y,8000,8,'E:\noise enhancing\wav\s&w_afterhf(10).wav');%0db带噪信号

frame = 256;    % Defining frame size

shift=128;
win=hamming(256);

%%%%%%     preemphasize     %%%%%%%%

signal(1)=y(1);
pre_u=0.8;
for j1 = 2:length(y),
 signal(j1) = y(j1)-pre_u*y(j1-1);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% for j1 = 1:length(y),
%  signal(j1) = y(j1);
% end;

hh = 0; 
head=0;
   
%%%%  estimation of noise  %%%%
   for k = 1 : 5,
        for m = 1 : frame,
          abc1(m) = signal(head+m);
        end;
        abc1=abc1.*win';
         
        head = head +frame;
        frame_temp(k,1:frame) = abs(fft(abc1));% FFT OF THE SIGNAL + NOISE FRAME BY FRAME
        frame_angle(k,1:frame) = angle(fft(abc1));% ANGLE OF FFT OF THE SIGNAL + NOISE FRAME BY FRAME
        ps_noise(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)));

    g=zeros(1,18);

    for i1=1:3,
	    g(1)=g(1)+ps_noise(k,i1);       
    end;
	Ec(k,1)=g(1)/(2+1);
    for i1=4:6,
		g(2)=g(2)+ps_noise(k,i1);   
    end;
    Ec(k,2)=g(2)/(2+1);
	for i1=7:10,
		g(3)=g(3)+ps_noise(k,i1);       
    end;
    Ec(k,3)=g(3)/(3+1);
	for i1=11:13,
		g(4)=g(4)+ps_noise(k,i1);    
    end;
    Ec(k,4)=g(4)/(2+1);
	for i1=14:16,
		g(5)=g(5)+ps_noise(k,i1);
    end;   
    Ec(k,5)=g(5)/(3+1);
	for i1=17:20,
		g(6)=g(6)+ps_noise(k,i1);
    end; 
    Ec(k,6)=g(6)/(3+1);
	for i1=21:25,
		g(7)=g(7)+ps_noise(k,i1);
    end;
	Ec(k,7)=g(7)/(4+1);
    for i1=26:29,
		g(8)=g(8)+ps_noise(k,i1);  
    end;
    Ec(k,8)=g(8)/(3+1);  
	for i1=30:35,
		g(9)=g(9)+ps_noise(k,i1);   
    end;
    Ec(k,9)=g(9)/(5+1);
	for i1=36:41,
		g(10)=g(10)+ps_noise(k,i1);    
    end;
    Ec(k,10)=g(10)/(5+1);
	for i1=42:47,
		g(11)=g(11)+ps_noise(k,i1);   
    end;
    Ec(k,11)=g(11)/(5+1);
	for i1=48:55,
		g(12)=g(12)+ps_noise(k,i1);    
    end;
    Ec(k,12)=g(12)/(7+1);
	for i1=56:64,
		g(13)=g(13)+ps_noise(k,i1);      
    end;
    Ec(k,13)=g(13)/(8+1);
	for i1=65:74,
		g(14)=g(14)+ps_noise(k,i1);   
    end; 
    Ec(k,14)=g(14)/(9+1);
	for i1=75:86,
		g(15)=g(15)+ps_noise(k,i1);   
    end;  
    Ec(k,15)=g(15)/(11+1);
	for i1=87:101,
		g(16)=g(16)+ps_noise(k,i1);  
    end;
    Ec(k,16)=g(16)/(14+1);
	for i1=102:118,
		g(17)=g(17)+ps_noise(k,i1);
    end;
    Ec(k,17)=g(17)/(16+1);
	for i1=119:128,
		g(18)=g(18)+ps_noise(k,i1);  
    end;
     Ec(k,18)=g(18)/(9+1);
% 	for i1=142:170,
% 		g(19)=g(19)+ps_noise(k,i1);   
%     end;Ec(k,19)=g(19)/(900+1);
% 	for i1=171:205,
% 		g(20)=g(20)+ps_noise(k,i1);  
%     end;
%     Ec(k,20)=g(20)/(1100+1);
% 	for i1=206:246,
% 		g(21)=g(21)+ps_noise(k,i1);   
%     end;  
%     Ec(k,21)=g(21)/(1300+1);
%     for i1=247:256,
% 		g(22)=g(22)+ps_noise(k,i1);    
%     end;
%     Ec(k,22)=g(22)/(1800+1);    
end

for i2=1:18,
   En(i2)=0;
   for k=1:5, 
       En(i2)=En(i2)+Ec(k,i2);      
   end
    En(i2)=En(i2)/5;
   if En(i2)<0.0625
       En(i2)=0.0625;
   end
end


head = 0; 
mm=1;
nn=1;
% START OF THE NOISE ELIMINATION THROUGH SPECTRAL SUBTRACTION BASED ON THE THRESHOLD SET

   for k = 1 : length(signal)/shift-1,
        for m = 1 : frame,
          abc1(m) = signal(head+m);  
        end;
        abc1=abc1.*win';
        
        head = head +shift;
        frame_temp(k,1:frame) = abs(fft(abc1));% FFT OF THE SIGNAL + NOISE FRAME BY FRAME

        frame_angle(k,1:frame) = angle(fft(abc1));% ANGLE OF FFT OF THE SIGNAL + NOISE FRAME BY FRAME
        
        ps_signal(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)));

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    g=zeros(1,18);
    sum_En(k)=0;
		
    for i1=1:3,
	    g(1)=g(1)+ps_signal(k,i1);       
    end;
	Ec(k,1)=g(1)/(2+1);
    for i1=4:6,
		g(2)=g(2)+ps_signal(k,i1);   
    end;
    Ec(k,2)=g(2)/(2+1);
	for i1=7:10,
		g(3)=g(3)+ps_signal(k,i1);       
    end;
    Ec(k,3)=g(3)/(3+1);
	for i1=11:13,
		g(4)=g(4)+ps_signal(k,i1);    
    end;
    Ec(k,4)=g(4)/(2+1);
	for i1=14:16,
		g(5)=g(5)+ps_signal(k,i1);
    end;   
    Ec(k,5)=g(5)/(3+1);
	for i1=17:20,
		g(6)=g(6)+ps_signal(k,i1);
    end; 
    Ec(k,6)=g(6)/(3+1);
	for i1=21:25,
		g(7)=g(7)+ps_signal(k,i1);
    end;
	Ec(k,7)=g(7)/(4+1);
    for i1=26:29,
		g(8)=g(8)+ps_signal(k,i1);  
    end;
    Ec(k,8)=g(8)/(3+1);  
	for i1=30:35,
		g(9)=g(9)+ps_signal(k,i1);   
    end;
    Ec(k,9)=g(9)/(5+1);
	for i1=36:41,
		g(10)=g(10)+ps_signal(k,i1);    
    end;
    Ec(k,10)=g(10)/(5+1);
	for i1=42:47,
		g(11)=g(11)+ps_signal(k,i1);   
    end;
    Ec(k,11)=g(11)/(5+1);
	for i1=48:55,
		g(12)=g(12)+ps_signal(k,i1);    
    end;
    Ec(k,12)=g(12)/(7+1);
	for i1=56:64,
		g(13)=g(13)+ps_signal(k,i1);      
    end;
    Ec(k,13)=g(13)/(8+1);
	for i1=65:74,
		g(14)=g(14)+ps_signal(k,i1);   
    end; 
    Ec(k,14)=g(14)/(9+1);
	for i1=75:86,
		g(15)=g(15)+ps_signal(k,i1);   
    end;  
    Ec(k,15)=g(15)/(11+1);
	for i1=87:101,
		g(16)=g(16)+ps_signal(k,i1);  
    end;
    Ec(k,16)=g(16)/(14+1);
	for i1=102:118,
		g(17)=g(17)+ps_signal(k,i1);
    end;
    Ec(k,17)=g(17)/(16+1);
	for i1=119:128,
		g(18)=g(18)+ps_signal(k,i1);  
    end;
     Ec(k,18)=g(18)/(9+1);
% 	for i1=142:170,
% 		g(19)=g(19)+frame_temp(k,i1);   
%     end;Ec(k,19)=g(19)/(900+1);
% 	for i1=171:205,
% 		g(20)=g(20)+ps_signal(k,i1);  
%     end;
%     Ec(k,20)=g(20)/(1100+1);
% 	for i1=206:246,
% 		g(21)=g(21)+ps_signal(k,i1);   
%     end;  
%     Ec(k,21)=g(21)/(1300+1);
%     for i1=247:256,
% 		g(22)=g(22)+ps_signal(k,i1);    
%     end;
%     Ec(k,22)=g(22)/(1800+1);
	
    for i2=1:18,
	  if k==1
        E_ch(k,i2)=Ec(k,i2);
      else  E_ch(k,i2)=0.45*E_ch(k-1,i2)+0.55*Ec(k,i2);
      end
      if E_ch(k,i2)<0.0625
         E_ch(k,i2)=0.0625;
      end
    end   
    
    for i2=1:18,
      sum_En(k)=sum_En(k)+En(i2);
    end
    
    gama_n(k)=max(-13,-10*log10(sum_En(k)));
    v_matrix=[2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 4 4 4 5 5 5 6 6 7 7 7 8 8 9 9 10 10 11 12 12 13 13 14 15 15 16 17 17 18 19 20 20 21 22 23 24 24 25 26 27
        28 28 29 30 31 32 33 34 35 36 37 37 38 39 40 41 42 43 44 45 46 47 48 49 50 50 50 50 50 50 50 50 50 50];
    v(k)=0;
    
    for i2=1:18,
        snr_ch(k,i2)=round(10*log10(E(k,i2)/En(i2))/0.375);
        delta(k,i2)=max(0,min(89,snr_ch(k,i2)));
        E_dB(k,i)=10*log10(Ec(k,18));
  %       gama_db(i2)=0.39*(delta(k,i2)-6)+gama_n(k);
 %       gama(i2)=min(1,power(10,gama_db(i2)/20));
    end
    
    E_dB_aver(k)=E_dB(k);
    
    for i2=1:18,
        v(k)=v(k)+v_matrix(delta(k,i2));
    end
    
    
    
    UPDATE_THLD=35;
    NOISE_FLOOR_DB=0;
    DEV_THLD=28;
    UPDATE_CNT_THLD=50;
    HYESTER_CNT_THLD=6;
    
    %% Normal update logic %%
    update_flag=1;        %%% false:1  true:0
    if ( v(k)<=UPDATE_THLD){
        update_flag=0;
        update_cnt=0;
       }
  
    %% Forced update logic %%
    else if ((E_tot(k)>NOISE_FLOOR_DB) and (DELTA_E<DEV_THLD))
         {
         update_cnt = update_cnt + 1;
         if(update_cnt=UPDATE_CNT_THLD)
             update_flag=0;
         end ;
         }
         end ;
     end
     %% "Hysteresis" logic to prevent long-term creeping of update_cnt %%
     if(update_cnt==last_update_cnt)
         hyster_cnt=hyster_cnt+1;
     else hyster_cnt=0;
     end ;
     last_update_cnt=update_cnt;
     if(hyster_cnt>HYESTER_CNT_THLD)
         update_cnt=0;
     end

         

    
    for i1=1:3,
	    frame1(k,i1)=frame_temp(k,i1)*gama(1);
    end;
	for i1=4:6,
        frame1(k,i1)=frame_temp(k,i1)*gama(2);
    end;
	for i1=7:10,
        frame1(k,i1)=frame_temp(k,i1)*gama(3);
    end;
	for i1=11:13,
		frame1(k,i1)=frame_temp(k,i1)*gama(4);
    end;
	for i1=14:16,
		frame1(k,i1)=frame_temp(k,i1)*gama(5);
    end;
	for i1=17:20,
		frame1(k,i1)=frame_temp(k,i1)*gama(6);
    end;
	for i1=21:25,
		frame1(k,i1)=frame_temp(k,i1)*gama(7);
    end;
	for i1=26:29,
		frame1(k,i1)=frame_temp(k,i1)*gama(8);
    end;
	for i1=30:35,
		frame1(k,i1)=frame_temp(k,i1)*gama(9);
    end;
	for i1=36:41,
		frame1(k,i1)=frame_temp(k,i1)*gama(10);
    end;
	for i1=42:47,
		frame1(k,i1)=frame_temp(k,i1)*gama(11);
    end;
	for i1=48:55,
		frame1(k,i1)=frame_temp(k,i1)*gama(12);
    end;
	for i1=56:64,
		frame1(k,i1)=frame_temp(k,i1)*gama(13);
    end;
	for i1=65:74,
		frame1(k,i1)=frame_temp(k,i1)*gama(14);
    end;
	for i1=75:86,
		frame1(k,i1)=frame_temp(k,i1)*gama(15);
    end;
	for i1=87:101,
		frame1(k,i1)=frame_temp(k,i1)*gama(16);
    end;
	for i1=102:118,
		frame1(k,i1)=frame_temp(k,i1)*gama(17);
    end;
	for i1=119:129,
		frame1(k,i1)=frame_temp(k,i1)*gama(18);
    end;
% 	for i1=142:170,
% 		frame1(k,i1)=frame_temp(k,i1)*gama(19);
%     end;
% 	for i1=171:205,
% 		frame1(k,i1)=frame_temp(k,i1)*gama(20);
%     end;
% 	for i1=206:246,
% 		frame1(k,i1)=frame_temp(k,i1)*gama(21);
%     end;    
%     for i1=247:256,
% 		frame1(k,i1)=frame_temp(k,i1)*gama(22);
%     end;
	
     for i1=130:256
     frame1(k,i1)=conj(frame1(k,258-i1));
     end
    
       frame1(k,1:frame) = frame1(k,1:frame).*(exp(i*frame_angle(k,1:frame)));
       frame2(k,1:frame)=ifft(frame1(k,1:frame));
       %signal(1,(((k-1)*frame)+1):(k*frame)) = frame2(k,1:frame); % Retriving back the signal(after spectral subtraction)
       if k==1 
             signal(1,1:shift)=frame2(k,1:shift);
       else   signal(1,((k-1)*shift+1):(k*shift))=(frame2(k,1:shift)+frame2(k-1,(shift+1):(shift*2)))/2;      
       end
   end;
   
%%%%%%%%%%  depreemphasize    %%%%%%%%%%%%%
signal2(1)=signal(1);
for j1 = 2:length(signal),
   signal2(j1) = signal(j1)+pre_u*signal2(j1-1);
end;
signal=signal2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  signal((length(y)-1000):(length(y)))=[];   %give up 4 frames in the end
  y1((length(y1)-1000):(length(y1)))=[];
%    y((length(y)-1000):(length(y)))=[];
% figure(3);
% plot(1:length(frame_ps),frame_ps,1:length(ps_final),ps_final);
figure(4);
signal=signal';
signal=signal/max(abs(signal));%归一化
plot(1:length(signal),signal);
%以下画语谱图
%map=(log10(1:0.1:10))';
%map=[map map map ];
%subplot(2,2,1)
%spgrambw( noise,8000);title('noise Specgram');colormap(1-map);
%subplot(2,2,2)
%spgrambw(y1,8000);title('sig Specgram');colormap(1-map);
%subplot(2,2,3)
%spgrambw(y,8000);title('mixsig Specgram');colormap(1-map);
%subplot(2,2,4)
%spgrambw(signal,8000);title('Proposed Algorithm Specgram');colormap(1-map);

%spgrambw(sig,8000);title('pure Specgram');colormap(1-map);
%[overall_snr2,seg_nr2]=snr(y1,signal)
% before_snr = 10*log10(sum(abs(y1).^2)/sum((abs(y1-y)).^2))
 overall_snr = 10*log10(sum(abs(y1).^2)/sum((abs(y1-signal)).^2))
wavwrite(signal,8000,8,'E:\noise enhancing\wav\mymasking_result(5).wav');