zm01xg_5_ed1.m

matlab仿真通过的降噪程序

clear all;
%  [y1,fs,bits]=wavread('F:\noise enhancing\wav\5.wav');
%  y1=y1/max(abs(y1));%语音信号归一化
%  wavwrite(y1,8000,8,'F:\noise enhancing\voise\5.wav');
%  figure(1);
%  plot(y1);
% 
%  [noise,fs1,bits1]=wavread('F:\noise enhancing\wav\5_noise.wav');
%  y=mixsig(y1,noise,10);% 混合
%  y=y/max(abs(y));%归一化
%  wavwrite(y,8000,8,'F:\noise enhancing\wav\mymasking_s&w(10).wav');
%  figure(2);
%  plot(y);

[y,fs,bits]=wavread('F:\noise enhancing\wav\radio_baby.wav');   %mys&w16bit.F:\noise enhancing\wav\wav  zhuiqiu_s&w(10) mys&w16bit
frame = 256;    % Defining frame size
shift=128;
win=hamming(256);

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

ps_noise=zeros(length(signal)/frame,frame);
frame_temp = zeros(length(signal)/frame,frame);

hh = 0; 
   for k = 1 : 5,
       for l = 1 : frame,
          b(l) = signal(hh+l);
      end;
        hh = hh + frame;
        frame_temp(k,1:frame) = abs(fft(b));                      %fft for the first 50 frames
        ps_noise(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)))/frame;
        %ps_noise(1,1:frame)= (sum(ps_noise(1:k,l))/20);          % Sum of the power spectral densities of samples within a frame

    end;%语音前五真早声能量
    %ps_noise=zeros(length(signal)/frame,frame);
    %ps_noise(1,1:frame)= sum(A)/20;            % setting the threshold for the noise(frame noise)
    ps_noise(1,1:frame)= (sum(ps_noise(1:k,1:frame))/5);
    frame_temp_initial(1:frame)=sum(frame_temp(k,1:frame))/k;

%
head = 0; 
mm=1;
nn=1;

da=0;xiao=0;equ=0;counter=0;

%%%%%% 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

        pmix_signal(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)))./frame; % power of mixed signal
        
        ps_temp(k,1:frame)=pmix_signal(k,1:frame);
        
        
%         c_n=0.995;   %%% rate of renewing noise power spectrum
%         c_s=0.999;    %%% rate of renewing signal power spectrum
%         c_beta=0.1;
%         
%         if k==1
%            ps_noise(k,1:frame)=c_n*ps_noise(1,1:frame)+(1-c_n)*pmix_signal(k,1:frame);
%            ps_signal(k,1:frame)=pmix_signal(k,1:frame);%%-c_beta*ps_noise(k,1:frame);
%         else
%            ps_noise(k,1:frame)=ps_noise(1,1:frame);
%            pmax(k)=max(pmix_signal(k,1:frame));
%            ps_signal(k,1:frame)=c_s*ps_signal(k-1,1:frame)+(1-c_s)*pmix_signal(k,1:frame)-c_beta*ps_noise(k,1:frame);
%            ps_signal(k,1:frame)=abs(ps_signal(k,1:frame));
%            % ps_signal(k,1:frame)=c_s*ps_signal(k-1,1:frame)+(1-c_s)*pmax(k)-c_beta*ps_noise(k,1:frame);
%         end
           
%%%%%%%%%%     add Cepstrum Distance   ( end detection ) %%%%%%%%%%%%%%%
       if k==1
             framenoise_temp(k,1:frame)=frame_temp_initial(1:frame);
       else 
             framenoise_temp(k,1:frame)=framenoise_temp(k-1,1:frame);
       end
       dd(k)=sum(ifft(log(abs(framenoise_temp(k,1:frame)))))/frame;  %c'噪声倒谱系数
       ddi=sum(ifft(log(abs(framenoise_temp(1,1:frame)))))/frame;
       d(k)=sum(pmix_signal(k,1:frame).*exp(-j*2*pi*(1:frame)*k/frame));%信号倒谱
       di=sum(pmix_signal(k,1:frame));
       
       if k>1 
           dd(k)=0.8*dd(k)+0.2*d(k-1);
       end
       
       cep_disp(k)=4.3429*sqrt((di-ddi).^2+2.*((d(k)-dd(k)).^2));%倒谱距离
       
       if k>1
           cep_disp(k)=0.8*cep_disp(k)+0.2*cep_disp(k-1);
       end
       
       c_n=0.995;   %%% rate of renewing noise power spectrum
       c_s=0.999;    %%% rate of renewing signal power spectrum
       c_beta=0.1; 
       
%        if cep_disp(k)>6  
%            counter=counter+1;
%        end
       
       if k==1
           ps_signal(k,1:frame)=pmix_signal(k,1:frame);
       else
         if cep_disp(k)<6  
           counter=counter+1;
           if counter==3 
               ps_noise(k,1:frame)=c_n*ps_noise(1,1:frame)+(1-c_n)*pmix_signal(k,1:frame);
               framenoise_temp(k,1:frame)=c_n*framenoise_temp(k-1,1:frame)+(1-c_n)*frame_temp(k,1:frame);
               counter=2;
           else ps_noise(k,1:frame)=ps_noise(k-1,1:frame);
           end

           
         else ps_noise(k,1:frame)=ps_noise(k-1,1:frame);
              counter=0;
             %   ps_signal(k,1:frame)=0;
         end
          ps_signal(k,1:frame)=c_s*ps_signal(k-1,1:frame)+(1-c_s)*pmix_signal(k,1:frame)-c_beta*ps_noise(k,1:frame);
       end

      ps_signal(k,1:frame)=abs(ps_signal(k,1:frame));
   
%         frame_ps(1,k) = (sum(ps_signal(k,1:frame)));%一阵内信号能量之和
%         frame_pn(1,k)=sum(ps_noise(k,1:frame));
%         %ps_final(1,k) = frame_ps(1,k) - threshold;
%         %aa=0.8;
%         ps_final(1,k) = frame_ps(1,k)- 0.8*frame_pn(1,k);
        
%      %短点检测
%         yeta(1,k)=sum(    ps_signal(k,1:frame)/(abs(ps_signal(k,1:frame)-ps_noise(k,1:frame))) );
%         %yeta(1,k)=frame_ps(1,k)/(frame_ps(1,k)- frame_pn(1,k));
%        if(yeta(1,k)>40),
%            ps_final(1,k) =0.1;
%        else 
%            if(yeta(1,k)>1)
%                ps_final(1,k)=(abs(frame_ps(1,k)-frame_pn(1,k)));
%            else
%            if ps_final(1,k)< 0,
%             ps_final(1,k) =0.1;
%         else
%              ps_final(1,k) = ps_final(1,k); 
%              
%          end;
%      end;
%        end;
       
%         aa=0.7;bb=2;
%         h(1,k)=power(ps_final(1,k)/(ps_final(1,k)+aa*frame_pn(1,k)),bb);
%         frame1(k,1:frame) = h(1,k).*(frame_temp(k,1:frame));             %%%%%%   Weiner滤波
%         
%         ps_signal(k,1:frame) = (frame1(k,1:frame).*conj(frame1(k,1:frame)))./frame;]
     
     
 %%%%%% calculate masking value   %%%%%%

	T=zeros(1,256);
    b=zeros(1,22);c=zeros(1,22);o=zeros(1,22);
    sf=zeros(22,22);
		
	for i1=1:3,
	    b(1)=b(1)+ps_signal(k,i1);
    end;
	for i1=4:6,
		b(2)=b(2)+ps_signal(k,i1);
    end;
	for i1=7:10,
		b(3)=b(3)+ps_signal(k,i1);
    end;
	for i1=11:13,
		b(4)=b(4)+ps_signal(k,i1);
    end;
	for i1=14:16,
		b(5)=b(5)+ps_signal(k,i1);
    end;
	for i1=17:20,
		b(6)=b(6)+ps_signal(k,i1);
    end;
	for i1=21:25,
		b(7)=b(7)+ps_signal(k,i1);
    end;
	for i1=26:29,
		b(8)=b(8)+ps_signal(k,i1);
    end;
	for i1=30:35,
		b(9)=b(9)+ps_signal(k,i1);
    end;
	for i1=36:41,
		b(10)=b(10)+ps_signal(k,i1);
    end;
	for i1=42:47,
		b(11)=b(11)+ps_signal(k,i1);
    end;
	for i1=48:55,
		b(12)=b(12)+ps_signal(k,i1);
    end;
	for i1=56:64,
		b(13)=b(13)+ps_signal(k,i1);
    end;
	for i1=65:74,
		b(14)=b(14)+ps_signal(k,i1);
    end;
	for i1=75:86,
		b(15)=b(15)+ps_signal(k,i1);
    end;
	for i1=87:101,
		b(16)=b(16)+ps_signal(k,i1);
    end;
	for i1=102:118,
		b(17)=b(17)+ps_signal(k,i1);
    end;
	for i1=119:141,
		b(18)=b(18)+ps_signal(k,i1);
    end;
	for i1=142:170,
		b(19)=b(19)+ps_signal(k,i1);
    end;
	for i1=171:205,
		b(20)=b(20)+ps_signal(k,i1);
    end;
	for i1=206:246,
		b(21)=b(21)+ps_signal(k,i1);
    end;    
    for i1=247:256,
		b(22)=b(22)+ps_signal(k,i1);
    end;
	
	
%	for i1=1:22,
%		for j1=1:22,
%		sf(i1,j1)=15.81+7.5*((i1-j1)+0.474)-17.5*sqrt(1+((i1-j1)+0.474)*((i1-j1)+0.474));    %$ (18)
%		sf(i1,j1)=power(10,sf(i1,j1)/(20));                                                  %$ transform from dB
%        end;
%    end;
    
    for i1=1:22,
    sf(i1)=15.81+7.5*(i1+0.474)-17.5*sqrt(1+(i1+0.474)*(i1+0.474));
    end;
		
%	for j1=1:22,
%		c(j1)=0;
%		for i1=1:22,
%			c(j1)=c(j1)+b(i1)*sf(i1,j1);                          %$ (19)
%        end;
%    end;
    
    cc_temp=conv2(b,sf);
    for i1=1:22		
		for j1=1:22
			c(i1)=c(i1)+cc_temp(i1,j1);