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

frame = 256;    % Defining frame size
shift=64;
win=hamming(256);

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

ps_noise=zeros(length(signal)/shift,frame);
frame_temp = zeros(length(signal)/shift,frame);
ps_temp= zeros(length(signal)/shift,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;
   end;
   ps_noise(1,1:frame)= (sum(ps_noise(1:5,1:frame))/5);

head = 0; 

%%%%%%%   Weiner Filtering  %%%%%%%%

   for k = 1 :( length(signal)/shift-3),
        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)))./frame;
        ps_noise(k,1:frame)=ps_noise(1,1:frame);
       
      frame_ps(1,k) = (sum(ps_signal(k,1:frame)));%一陣內信號能量之和
      frame_pn(1,k) = sum(ps_noise(k,1:frame));
%       if k>3
%           ps_final(1,k) = (frame_ps(1,k-2)+frame_ps(1,k-1)+frame_ps(1,k))/3- 0.8*frame_pn(1,k); 
%       else
%           ps_final(1,k) = frame_ps(1,k)- frame_pn(1,k);  
%       end
       
      ps_final(1,k) = frame_ps(1,k)- 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 
       
%        ps_temp(k,1:frame)=ps_signal(k,1:frame)-ps_noise(k,1:frame);
%        if ps_temp(k,1:frame)<0
%            ps_temp(k,1:frame)=0;
%        end
%       frame_temp(k,1:frame)=sqrt( ps_temp(k,1:frame));
       aa=0.7;bb=2;
       h(1,k)=power(ps_final(1,k)/(ps_final(1,k)+aa*frame_pn(1,k)),bb);
       frame_temp(k,1:frame) = h(1,k).*(frame_temp(k,1:frame));%為那濾波
 
       
       frame_temp(k,1:frame) = frame_temp(k,1:frame).*(exp(i*frame_angle(k,1:frame)));
       frame2(k,1:frame)=ifft(frame_temp(k,1:frame));
       if k==1 
            signal(1,1:shift)=frame2(k,1:shift);
        else if k==2
              signal(1,(shift+1):(2*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 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((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);

overall_snr = 10*log10(sum(abs(y1).^2)/sum((abs(y1-signal)).^2))
wavwrite(signal,8000,8,'F:\noise enhancing\wav\SS\result(10db).wav');