svd_sine.m

用奇异值分解来消除正弦信号的噪声

   clear all ; close all;
   M=300;
   s_ori=sin(2*pi*(1:M)/40)'+2*cos(2*pi*(1:M)/30)'+cos(2*pi*(1:M)/60)';
   sig=awgn(s_ori,15,'measured');
   xds2=sig;

%%%%%%%%%%%%%%%以下内容是我自己所加，目的是对ca2=xds2进行奇异值分解去噪
 N=100;    ca2=xds2;
     ca2_new=zeros(size(ca2));
for  i=1:(M/100)
     w=zeros(N,1);
     w=ca2(100*(i-1)+1:100*i);
     ww=zeros(10,10);
for  num=1:10
     ww(num,:)=w(1+10*(num-1):10*num)';
end;
     [u,s,v]=svd(ww);  ss=s;
     total=TRACE(s);
     sigma_total=0;
for  p=1:10
    if sqrt(sigma_total/total)>=0.85
        break;
    else  sigma_total=sigma_total+s(p,p);
    end;
end;
       p_end=p;
       
for  p=p_end:10
     s(p,p)=0;
end;
     s_new=s;
     ww_new=u*s_new*v';
     w_new=zeros(N,1);
     
for  num=1:10
     w_new(10*(num-1)+1:10*num)=ww_new(num,:)';
end;
     ca2_new(100*(i-1)+1:100*i)=w_new;
end;
   
    xds2_new=ca2_new;
    s_new=xds2_new;
    %%%%%   所增加的内容完毕
    %%%%  观看矩阵的奇异值
        ss_sigular=zeros(10,1);
    for i=1:10
        ss_sigular(i)=ss(i,i);
    end;
        s_sigular=zeros(10,1);
    for i=1:10
        s_sigular(i)=s(i,i);
    end;
    figure;
    subplot(211);
    stem(ss_sigular); title('原始信号奇异值'); grid on;
    subplot(212);
    stem(s_sigular,'r');  title('消除噪声后的奇异值');    grid on;
    
     figure;
  subplot(311);
  plot(s_ori); axis tight; grid on;
  subplot(312);
  plot(sig);axis tight;   grid on;
  subplot(313);
  plot(s_new,'r');axis tight; grid on;
  -1*10*log10((s_new-s_ori)'*(s_new-s_ori)/(s_ori'*s_ori))