%%% @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ %%%
 
%     This algorithm is used to enhance the noisy speech especially poluted by broadband white noise.
%     
%     Last modified  April 2004
%     
%     Southeast University

%%% @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ %%%

clear all;
[y1,fs,bits]=wavread('E:\noise enhancing\wav\5.wav');      % input clean speech
y1=y1/max(abs(y1));                                                 % normalize
figure(1);
plot(y1);                                                           % depict clean speech

[noise,fs1,bits1]=wavread('E:\noise enhancing\wav\5_noise.wav');      % input noise
y=mixsig(y1,noise,0);                                                          % mix clean speech with noise according to definite SNR
y=y/max(abs(y));                                                               % renormalize
%wavwrite(y,8000,8,'F:\noise enhancing\wav\mymasking_s&w(0).wav');              % retrieve noisy speech
figure(2);
plot(y);

%% [y,fs,bits]=wavread('F:\noise enhancing\noise\denoise\stc-launch[1].wav'); % or input the noisy speech
frame = 256;              % Defining frame size
shift=128;                % Defining frame shift
win=hamming(256);

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

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

%%%%  estimation of noise energy using the first five frames  %%%%
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 5 frames
        ps_noise(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)))/frame;   % power
        %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(1,1:frame)= (sum(ps_noise(1:k,1:frame))/5);          % average power spectral densities of noises within a frame



%%%%% START OF THE NOISE ELIMINATION  %%%%%%
head = 0; 
   for k = 1 :( length(signal)/shift-1),
        for m = 1 : frame,
          abc1(m) = signal(head+m);          
        end;
        
        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);
                  
%%%%%%%%%%%      end detection    %%%%%%%%%%%%%%%

       c_n=0.995;    %%% rate of renewing noise power spectrum
       c_s=0.999;    %%% rate of renewing signal power spectrum
       c_beta=0.1;  
       counter=0;    %%% the number of consecutive noise frames till current frame

     if k==1                                               % first frame
            ps_signal(k,1:frame)=pmix_signal(k,1:frame);
            counter=1;
     else  
       yeta(1,k)=sum(    pmix_signal(k,1:frame)/(abs(pmix_signal(k,1:frame)-ps_noise(k-1,1:frame))) );
       if(yeta(1,k)>1.1)                     %if yeta is larger than threshold,current frame is signed as noise frame
           counter=counter+1;
           if counter==3                     %if three consecutive frames are signed as noise frame,renew the estimation of noise energy
               ps_noise(k,1:frame)=c_n*ps_noise(1,1:frame)+(1-c_n)*pmix_signal(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;             
       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);  %%% estimation of power spectrum of signal
     end
       ps_signal(k,1:frame)=abs(ps_signal(k,1:frame));
        
     
%%%%%%%%%   calculate masking value   %%%%%%%%%

	T=zeros(1,128);
    b=zeros(1,18);c=zeros(1,18);o=zeros(1,18);
    sf=zeros(18,18);

 %%%  calculate power spectrum of every bark band    
    
	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:128,
		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;
	
%%%  calculate the spread function
    
    for i1=1:18,
    sf(i1)=15.81+7.5*(i1+0.474)-17.5*sqrt(1+(i1+0.474)*(i1+0.474));
    end;

%%%  apply the spread function to the critical band spectrum

%	 for j1=1:22,
%		c(j1)=0;
%		for i1=1:22,
%			c(j1)=c(j1)+b(i1)*sf(i1,j1);                         
%        end;
%    end;                                % the function of this section is same as the following function(conv2)
    
    cc_temp=conv2(b,sf);                 % convolution 
    for i1=1:18		
		for j1=1:18
			c(i1)=c(i1)+cc_temp(i1,j1);
        end;
    end;

%%% calculate the spread masking threshold

	temp_value(k)=0.0;
	for i1=1:18,
		temp_value(k)=temp_value(k)+b(i1);
    end;
	ua=temp_value(k)/128.0;               % calculate the arithmetic average

	temp_value(k)=0;
	for i1=1:128, 
		temp_value(k)=temp_value(k)+log10(ps_signal(k,i1));
    end;
	temp_value(k)=temp_value(k)/128.0;
	uj=power(10,temp_value(k));           % calculate the geomitry average
    
	 sfm=-10*log10(uj/ua);                % spectrual flatness measurement

	u=min(sfm/(-60),1);                   % tonality of signal

  for i1=1:18,
	  O(i1)=u*(14.5+i1)+(1-u)*5.5;            % offset
      T(i1)=power(10,log10(c(i1))-O(i1)/10);  % auditory masking threshold
	  c(i1)=T(i1);                            % T[i] temperorily stored in c[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:128,
		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;
	
		
%%%  calculate the absolute threshold
	mm=0.0;
	for i1=1:129,
		f(i1)=mm;
		mm=mm+8/129;
    end;
	
	f(1)=f(2);
	for i1=1:129,
         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:129,
		T(i1)=max(T(i1),f(i1));
    end;
    
    
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       denoise  algorithm        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

landa=1.5;          %%%  make up for estimate error of a_l & a_h  
c_a=0.7;            %%%  coefficient of a_1 & a_2
a(k,1:(frame/2))=(c_a*a1(k,1:(frame/2))+(1-c_a)*a2(k,1:(frame/2)))*landa; 

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

for i1=1:129,
  frame1(k,i1)=sqrt(e_ps_signal(k,i1));
end
%%%  the later 128 points
   for i1=130:256
     frame1(k,i1)=conj(frame1(k,258-i1));   
   end

%%%  add phase to the signal
%frame1(k,1:(frame)) = sqrt(e_ps_signal(k,1:(frame))).*(exp(i*frame_angle(k,1:(frame))));
frame1(k,1:(frame)) = frame1(k,1:(frame)).*(exp(i*frame_angle(k,1:(frame))));

%%%  back to time domain
frame2(k,1:frame)=ifft(frame1(k,1:frame));

%%%  Retriving back the signal
   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     
 

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);
signal=signal';
signal=signal/max(abs(signal));     % normalize
plot(1:length(signal),signal);

overall_snr = 10*log10(sum(abs(y1).^2)/sum((abs(y1-signal)).^2))   %%%  estimate the SNR after processing

wavwrite(signal,8000,8,'E:\noise enhancing\wav\result1(0).wav');