www.pudn.com > Svd__Sine.zip > 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))