end_fenxing.m

matlab仿真通过的降噪程序

clear all;
[y1,fs,bits]=wavread('f:\voise\3.wav');
y1=y1/max(abs(y1));%归一化
wavwrite(y1,8000,8,'f:\wav\3.wav');
figure(1);
plot(y1);

[noise,fs1,bits1]=wavread('f:\voise\3_noise.wav');
y=mixsig(y1,noise,-5);
y=y/max(abs(y));%归一化
wavwrite(y,8000,8,'f:\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);
for j1 = 1:length(y),
 signal(j1) = y(j1);
 signal_ori(j1)=y1(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);


error2=0;
error=0;    %number of points that is wrongly detected 
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);
          abc_ori(m)=signal_ori(head+m);
          %abc1=abc1.*win';
        end;
        
        m=1:(frame-1);
        dd1(k)=sum(abs(abc1(m)-abc1(m+1)));
        for m=1:(frame/2-1);
            maxnum(m)=max(abc1((m*2-1):(m*2+1)));
            minnum(m)=min(abc1((m*2-1):(m*2+1)));
        end
        m=1:(frame/2-1);
        dd2(k)=sum(maxnum(m)-minnum(m));
       
        df(k)=1+log2(dd1(k))-log2(dd2(k));                %分形
        
       % ori_value(k)=sum(abc_ori(m));
    %    if ori_value(k)>0
    %        ori(k)=0;
    %      else ori(k)=1;
    %      end
        
        if df(k)>1.42
            dec(k)=0;
        else dec(k)=1;
        end
        
        ori_value(k)=sum(abc_ori(m));
        if ori_value(k)==0
            ori(k)=0;
        else ori(k)=1;
        end
        
%        if dec(k)==ori(k)
%       else error=error+1;
%        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)=0.99*ps_noise(k-1,1:frame)+0.01*ps_signal(k,1:frame);
           ps_signal(k,1:frame)=0.2*ps_signal(k-1,1:frame)+0.8*ps_signal(k,1:frame);
        end
        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.75*frame_pn(1,k);
        
        if df(k)>1.36
            ps_final(1,k) = 0.01;
        else if ps_final(1,k)<0
                ps_final(1,k)=0.01;
            else  ps_final(1,k) = frame_ps(1,k)- 0.8*frame_pn(1,k);
            end
        end
      
    if dec(k)==ori(k)
    else error=error+1;
    end
        
    if dec(k)==0 & ori(k)==1
       error2=error2+1;
    end
   end
   figure(3);
   stem(ori);
   figure(4);
   stem(dec);