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.wav');%0db帶噪信號
figure(2);
plot(y);
%clear all;
%[y1,fs,bits]=wavread('8k8bit.wav'); % Actual Signal
%[noise,fs1,bits1]=wavread('8k8bit_noise.wav');
%y=mixsig(y1,noise, 10);
%wavwrite(y,8000,8,'f:\a\mymasking_s&w.wav');
%figure(1);
%plot(y);
frame = 256;    % Defining frame size
%for k = 1:2560, % Loop for first 50 frames(0.5 seconds) of noise
%  y_temp(k) = y(k);
%end;
shift=128;
win=hamming(256);
%a=length(y),a=18800
for j1 = 1:length(y),
 signal(j1) = y(j1);
end;

ps_noise=zeros(length(signal)/frame,frame);
frame_temp = zeros(length(signal)/frame,frame);
%length1 = length(y_temp);   % length of the noise samples(first 4000 samples)
%ps_noise = zeros(length1/frame,frame);
%frame_temp = zeros(length1/frame,frame);
%hh = 0; 
 %  for k = 1 : length1/frame,
 %       for l = 1 : frame,
%          b(l) = y_temp(hh+l);
%        end;
%        hh = hh + frame;
%        frame_temp(k,1:frame) = 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;
%        A(1,k) = (sum(ps_noise(k,1:frame)));          % Sum of the power spectral densities of samples within a frame

%  end;
%threshold = sum(A)*frame/length1;            % setting the threshold for the noise(frame noise)

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);



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

   for k = 1 : length(signal)/frame,
        for m = 1 : frame,
          abc1(m) = signal(head+m);
          %abc1=abc1.*win';
        end;
        
        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_signal(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)))./frame;
        
        ps_temp=zeros(1,frame);
        ps_temp(1,1:frame)=ps_signal(k,1:frame);
        
        
        if k==1
            ps_noise(k,1:frame)=0.98*ps_noise(1,1:frame)+0.02*ps_signal(k,1:frame);
           % ps_signal(k,1:frame)=0.98*ps_signal(k,1:frame)+0.02*ps_signal(k,1:frame);
        else
           ps_noise(k,1:frame)=ps_noise(1,1:frame);%0.99*ps_noise(1,1:frame)+0.01*ps_signal(k,1:frame);
           ps_signal(k,1:frame)=ps_signal(k,1:frame);%0.99*ps_signal(k,1:frame)+0.01*ps_noise(1,1:frame);%+0.01*ps_signal(k-1,1:frame);
        end
        frame_ps(1,k) = (sum(ps_signal(k,1:frame)));%一陣內(nèi)信號能量之和
        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;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


	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 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)=sum(cc_temp(i1,j1));
        end;
    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:1010,
		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;
        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;
        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));%為那濾波
        
        for k1=1:256,
			%ps_noise[k]=0.98*ps_noise[k]+0.02*ps_signal[k];
			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=threshold/256/abs(frame1(k,k1));
            
            temppow=abs(ps_noise(k,k1))/abs(frame1(k,k1));
			if(power(temppow,2)<(1.0/(alfa(k1)+beita(k1))))
				%if(1-alfa[k]*pow(temppow,2)>0)
				frame1( k,k1)=sqrt(power(1-alfa(k1)*power(temppow,2),0.5))*frame1(k,k1);
           else
				frame1(k,k1)=sqrt(power(beita(k1)*power(temppow,2),0.5))*frame1(k,k1);
            end;
        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)
        
   end;
 signal((length(y)-1000):(length(y)))=[];   %give up 4 frames in the end
 y1((length(y1)-1000):(length(y1)))=[];
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)
overall_snr = 10*log10(sum(abs(y1).^2)/sum((abs(y1-signal)).^2))
wavwrite(signal,8000,8,'F:\noise enhancing\wav\mymasking21(10db).wav');