zm02.m

matlab仿真通过的降噪程序

    end;
	
	temp_value=0.0;
	for i1=1:22,
		temp_value=temp_value+b(i1);
    end;
	ua=temp_value/256.0;

	temp_value=0;
	for i1=1:256, 
		temp_value=temp_value+log10(ps_signal(k,i1));
    end;
	temp_value=temp_value/256.0;
	uj=power(10,temp_value);
    
	 sfm=-10*log10(uj/ua);                %$ (20)

	u=min(sfm/(-60),1);                   %$ (21)

  for i1=1:22,
	  O(i1)=u*(14.5+i)+(1-u)*5.5;
      T(i1)=power(10,log10(c(i1))-O(i1)/10);
	  c(i1)=T(i1);%c[i]暂存T[i]
  end;
    for i1=1:3,
	    T(i1)=c(1);
    end;
	for i1=4:6,
		T(i1)=c(2);
    end;
	for i1=7:10,
		T(i1)=c(3);
    end;
	for i1=11:13,
		T(i1)=c(4);
    end;
	for i1=14:16,
		T(i1)=c(5);
    end;
	for i1=17:20,
		T(i1)=c(6);
    end;
	for i1=21:25,
		T(i1)=c(7);
    end;
	for i1=26:29,
		T(i1)=c(8);
    end;
	for i1=30:35,
		T(i1)=c(9);
    end;
	for i1=36:41,
		T(i1)=c(10);
    end;
	for i1=42:47,
		T(i1)=c(11);
    end;
	for i1=48:55,
		T(i1)=c(12);
    end;
	for i1=56:64,
		T(i1)=c(13);
    end;
	for i1=65:74,
		T(i1)=c(14);
    end;
	for i1=75:86,
		T(i1)=c(15);
    end;
	for i1=87:101,
		T(i1)=c(16);
    end;
	for i1=102:118,
		T(i1)=c(17);
    end;
	for i1=119:141,
		T(i1)=c(18);
    end;
	for i1=142:170,
		T(i1)=c(19);
    end;
	for i1=171:205,
		T(i1)=c(20);
    end;
	for i1=206:246
		T(i1)=c(21);
    end;
	for i1=247:256,
		T(i1)=c(22);
    end;
	
		
	%////计算绝对听阈////
	mm=0.0;
	for i1=1:256,
		f(i1)=mm;
		mm=mm+8/256;
    end;
	
	f(1)=f(2);
	for i1=1:256,
		 %f(i1)=(3.64*power(f(i1),-0.8)-6.5*exp(-0.6*power((f(i1)-3.3),2))+0.001*power(f(i1),4));
         f(i1)=3.64*(f(i1).^(-0.8))- 6.5*exp(-0.6*((f(i1)-3.3).^2))+0.001*(f(i1).^4);
     end;
	
	for i1=1:256,
		T(i1)=max(T(i1),f(i1));
    end;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       original  algorithm            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         tmax=0;
%         tmin=0;
% 	for i1=1:256,
% 		if(T(i1)>tmax)
% 			tmax=T(i1);
%         end;
% 		if(T(i1)<tmin)
% 			tmin=T(i1);
%         end;
%     end;
         			
        %frame_ps(1,k) = (sum(ps_signal(k,1:frame)));% Sum of the power spectral densities of samples within a frame
        %ps_final(1,k) = frame_ps(1,k) - threshold;% Elimination of noise from the corrupted signal

        %h(1,k)=abs(ps_final(1,k)/frame_ps(1,k));
        %aa=0.6;bb=2;
     
      
        
%         for k1=1:256,
% 			alfa(k1)=(tmax*6-T(k1)*6+1*T(k1)-1*tmin)/(tmax-tmin);
% 			beita(k1)=(tmax*0.02-T(k1)*0.02+0*T(k1)-0*tmin)/(tmax-tmin);      
%             temppow=abs(ps_noise(k,k1))/abs(ps_temp(k,k1));
%             
% 			if(power(temppow,2)<(1.0/(alfa(k1)+beita(k1))))
% 				%if(1-alfa[k]*pow(temppow,2)>0)
% 				frame1( k,k1)=power((1-alfa(k1)*power(temppow,2)),0.5)*ps_temp(k,k1);
%             else
% 				frame1(k,k1)=power((beita(k1)*power(temppow,2)),0.5)*ps_temp(k,k1);
%             end;
%         end;
%  end;       
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       new  algorithm        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% a1(k,1:frame)=power((ps_signal(k,1:frame)+ps_noise(k,1:frame)),2)./(power((sqrt(ps_signal(k,1:frame))+sqrt(T(1:frame))),2))-(ps_signal(k,1:frame)+ps_noise(k,1:frame));
% a2(k,1:frame)=power((ps_signal(k,1:frame)+ps_noise(k,1:frame)),2)./(power((sqrt(ps_signal(k,1:frame))-sqrt(T(1:frame))),2))-(ps_signal(k,1:frame)+ps_noise(k,1:frame));

a1(k,1:frame)=power((ps_signal(k,1:frame)+ps_noise(k,1:frame)),2)./(ps_signal(k,1:frame))-(ps_signal(k,1:frame)+ps_noise(k,1:frame));
a2(k,1:frame)=power((ps_signal(k,1:frame)+ps_noise(k,1:frame)),2)./(T(1:frame))-(ps_signal(k,1:frame)+ps_noise(k,1:frame));

if a1(k,1:frame)>a2(k,1:frame)
     a(k,1:frame)=a1(k,1:frame);
else a(k,1:frame)=a2(k,1:frame);
end

landa=5;      %%%  make up for estimate error of a_l & a_h  
c_a=0.4;            %%%  coefficient of a_l & a_h
%a(k,1:frame)=(c_a*a_l(k,1:frame)+(1-c_a)*a_h(k,1:frame))*landa; 
%a(k,1:frame)=(ps_noise(k,1:frame)+T(1:frame)).*(ps_noise(k,1:frame)/T(1:frame));

e_ps_signal(k,1:frame)=power(pmix_signal(k,1:frame),2)./(pmix_signal(k,1:frame)+a(k,1:frame));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         Weiner Filter      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    hh = 0; 
%    for k = 1 : 5,   
%            ps_noise(k,1:frame) = (frame1(k,1:frame).*conj(frame1(k,1:frame)))/frame;
%    end;
%    %语音前五真早声能量
%    ps_noise2(1,1:frame)= (sum(ps_noise(1:5,1:frame))/5);
%   % frame1_initial(1:frame)=sum(frame1(1:5,1:frame))/5;
%    
%    for k = 1 : (length(signal)/shift-3),
%       ps_noise(k,1:frame)=ps_noise2(1,1:frame);
%       ps_signal(k,1:frame) = (frame1(k,1:frame).*conj(frame1(k,1:frame)))./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)-0* frame_pn(1,k);      
%        
%       %ps_final(k,1:frame) = ps_signal(k,1:frame)-ps_noise(k,1:frame);
%       %if ps_final(k,1:frame)<0
%       %    ps_final(k,1:frame)=zeros(1,frame);
%       %end
%       
%       %frame1(k,1:frame)=sqrt(ps_final(k,1:frame));
%       
%       if ps_final(1,k)<0
%           ps_final(1,k)=0;
%       end 
%       
%        aa=1;bb=0.5;
%        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).*(frame1(k,1:frame));%为那滤波
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

       frame1(k,1:frame) = sqrt(e_ps_signal(k,1:frame)).*(exp(i*frame_angle(k,1:frame)));
       frame2(k,1:frame)=ifft(frame1(k,1:frame));
        if k==1 
             signal(1,1:shift)=frame2(k,1:shift);
             %        else if k==2
        else   signal(1,((k-1)*shift+1):(k*shift))=(frame2(k,1:shift)+frame2(k-1,(shift+1):(shift*2)))/2;
%         else if  k==3
%               signal(1,(shift*2+1):(3*shift))=(frame2(k,1:shift)+frame2(k-1,(shift+1):(shift*2))+frame2(k-2,(shift*2+1):(shift*3)))/3;
%               %else if k==length(signal)/shift-3
%             %signal(1,(k*shift+1):((k+1)*shift))=frame2(k,(3*shift+1):frame);
%         else signal(1,(((k-1)*shift+1):(k*shift)))=(frame2(k,1:shift)+frame2(k-1,(shift+1):(shift*2))+frame2(k-2,(shift*2+1):(shift*3))+frame2(k-3,(shift*3+1):frame))/4;
%             %end;
%             end;
%            end;
         end;
       
       %signal(1,(((k-1)*frame)+1):(k*frame)) = frame2(k,1:frame); % Retriving back the signal(after spectral subtraction)
    end     
 

 signal((length(y)-1000):(length(y)))=[];   %give up 4 frames in the end
  signal(1:1280)=[];
  y1((length(y1)-1000):(length(y1)))=[];
 y1(1:1280)=[];
%figure(3);
%plot(1:length(frame_ps),frame_ps,1:length(ps_signal),ps_signal);
figure(4);
signal=signal';
signal=signal/max(abs(signal));%归一化
plot(1:length(signal),signal);

overall_snr = 10*log10(sum(abs(y1).^2)/sum((abs(y1-signal)).^2))
%overall_snr2 = 10*log10(sum(abs(y1).^2)/sum((abs(y1-signal*1.2)).^2))

wavwrite(signal,8000,8,'F:\noise enhancing\wav\mymasking21(5_odb).wav');

% for i=33:length(signal)
%     signal3(i)=signal(i)+signal(i-32);
% end
% 
% wavwrite(signal3,8000,8,'F:\noise enhancing\wav\mymasking21(5db).wav');

% for i=33:length(signal)
%     signal4(i)=signal(i)+signal(i-16);
% end
% wavwrite(signal4,8000,8,'F:\noise enhancing\wav\mymasking21(-5db).wav');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     以下画语谱图     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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)