svd_shocksin.m

自编的奇异值分解消噪程序

clc
clear all
close all

Ts=0.001;
num_t=500;
t=0:Ts:Ts*num_t;

n1=randn(numel(t),1);
n1=(n1-mean(n1));
vary1=var(n1);
n1=n1/sqrt(vary1);

z1=10*sin(2*pi*10*t);
z2=[zeros(1,90) cos(2*pi*100*t(1,91:130))  zeros(1,numel(t)-130)];

xx=z1+z2;
x=xx+max(z2)*n1.'*0;

%z1=2*sin(2*pi*50*t);
%z2=5*sin(2*pi*75*t);
%x=z1+z2+0*rand(1,numel(t));

%x=1*sin(100*pi*t)+0.3*sin(250*pi*t)+0.3*sin(500*pi*t);  % compare with SVD
%x=1*sin(100*pi*t)+5*sin(150*pi*t);
%x=10*cos(2*pi*100*t)+cos(2*pi*250*t);


%x=10*cos(2*pi*10*t)+[zeros(1,49) cos(2*pi*600*t(1,50:60))  zeros(1,numel(t)-60)];

%x=10*cos(2*pi*10*t)+cos(2*pi*60*t)+0.5*randn(1,numel(t));

%x=10*cos(2*pi*10*t)+[zeros(1,49) 10*cos(2*pi*600*t(1,50:60))  zeros(1,numel(t)-60)]+randn(1,numel(t));


num_x=numel(x);
num_row=round(numel(x)*0.9);

num_col=num_x-num_row+1;
A=hankel( x(1:num_row),x(num_row:num_x) );

[U E V]=svd(A);

%--------------------------------------------------
sv=diag(E); % all single value
figure(2)
semilogy(1:numel(sv),sv,'k-o')
xlabel('Number of singular values');
ylabel('Singular values');

%valve_e=0.8;
%ii=1;
%while sum( sv(1:ii) )/sum(sv)<=valve_e
%    ii=ii+1;
%end
%*************************************************************************
%rk=ii;
rk=2
%*************************************************************************
%SNR=sum( sv(1:rk) )/sum(sv(rk+1:rank(A))); % signal to noise ratio
%--------------------------------------------------
A_s=U(1:num_row,1:rk) * E(1:rk,1:rk) * V(1:num_col,1:rk).';

kk=1;
for ii=1:num_row
    x_s(kk)=A_s(ii,1);
    kk=kk+1;
end

for jj=2:num_col
    x_s(kk)=A_s(num_row,jj);
    kk=kk+1;
end

figure(3)
%plot(t,x,'k',t(1:numel(x_s)),x_s,'r')
plot(t,x,'k',t,x_s,'r')
%-----------------------------------------------------
%*************************************************************************

Ur=U;
Vr=V;

bb=1;
dd=1;
for rr=1:rk
    
row_x2=num_row;    
        
sv4=1;
sv2=1;

%------------------------------------------------------------------
while sv4>=sv2
    
x2=Ur(1:row_x2,rr);
num_x2=numel(x2);
row_x2=round(numel(x2)*0.9);
col_x2=num_x2-row_x2+1;

B=hankel( x2(1:row_x2),x2(row_x2:num_x2) );

clear U1 E1 V1
[U1 E1 V1]=svd(B);

sv_x2=diag(E1); % all single value

figure(4)
plot(sv_x2,'k-o')

sv2=sum(sv_x2(1:2));
sv4=sum(sv_x2(3:4));
rk_x2=2;

B_s=U1(1:row_x2,1:rk_x2) * E1(1:rk_x2,1:rk_x2) * V1(1:col_x2,1:rk_x2).';

clear x2_s;

kk=1;
for ii=1:row_x2
    x2_s(kk)=B_s(ii,1);
    kk=kk+1;
end

for jj=2:col_x2
    x2_s(kk)=B_s(row_x2,jj);
    kk=kk+1;
end


x2_s=x2_s.';
Ur(1:row_x2,rr)=x2_s(1:row_x2,1);

end
%------------------------------------------------------------------

row_min(bb)=row_x2;
bb=bb+1;

Ut(1:row_x2,dd)=x2_s(1:row_x2,1);
dd=dd+1;

end

%*************************************************************************

row_x2=min(row_min);
x_sum=0;
cc=1;
for qq=1:rr
    
    Es=E(qq,qq);
    Vs=V(:,qq);

    %for vv=1:bb
    Us=Ut(1:row_x2,cc);

    
    C_s=Us * Es * Vs.'; 

    pp=1;
    for ii=1:row_x2
    x_u(pp,cc)=C_s(ii,1);
    pp=pp+1;
    end
    for jj=2:num_col
    x_u(pp,cc)=C_s(row_x2,jj);
    pp=pp+1;
    end

    x_sum=x_sum+x_u(:,qq);
    
    cc=cc+1;
    
         
end

%--------------------------------------------------------------------------
% determine correlation
x_d=x_u(1:row_x2,:);
[row_x col_x]=size(x_d);
cx_valve=0.4;

ii=1;
jj=1;
kk=1;
tt=1;

x_cx(1:row_x,1)=0;
while ii<=col_x
    while jj<=col_x
        cx=mscohere(x_d(:,ii),x_d(:,jj));
        if mean(cx)>cx_valve & jj>=ii

            sign_jj(tt)=jj;tt=tt+1;
            [r_jj c_jj v_jj]=find(sign_jj==ii);
            if numel(c_jj)>1
               break;
            end
           
            x_cx(:,kk)=x_cx(:,kk)+x_d(:,jj);
            
        end
        jj=jj+1;
    end

    kk=kk+1;
      
    ii=ii+1;
    
    jj=1;
    x_cx(1:row_x,kk)=0;
end

[row_xcx col_xcx]=size(x_cx);
kk=1;
for ii=1:col_xcx
    if max(abs(x_cx(:,ii)))~=0
    x_last(:,kk)=x_cx(:,ii);
    kk=kk+1;
    end
end

%--------------------------------------------------------------------------
% Reconstruct data which has been truncated

[row_last col_last]=size(x_last);

for jj=1:col_last

% maxmium
max_ii=find(diff(sign(diff(x_last(:,jj))))==-2)+1;
max_val=x_last(max_ii,jj);

% minimum
min_ii=find(diff(sign(diff(x_last(:,jj))))==2)+1;
min_val=x_last(min_ii,jj);

% amplitude for the siganl
if numel(max_ii)==0 | numel(min_ii)==0
   disp('Length of the data is not enough! Results will be not right. Please add more data. Any button to continue.');pause
else
amp_val=(sum(abs(max_val))/numel(max_val)+sum(abs(min_val))/numel(min_val))/2;
end

% frequency
T0=(max_ii(numel(max_ii))-max_ii(1))*Ts/(numel(max_ii)-1);
f0=1/T0;

% harmonic signal

t0=t(max_ii(1));
x_rec1=amp_val*cos(2*pi*f0*(t-t0));
x_rec2(:,jj)=x_rec1(1:numel(t)).';

end


%--------------------------------------------------------------------------
%x_remain(:,1)=x(1:row_x).'-sum(x_last.').';
x_remain(:,1)=x.'-x_rec2;

%--------------------------------------------------------------------------


figure(10)
%plot(t(1:row_x2),x_last(:,1),'b',t(1:row_x2),z1(1:row_x2),'r:')
plot(t,x_rec2(:,1),'k',t,z1,'r:')
xlabel('Time/s');
ylabel('x_1');
legend('Decomposed','Original');
%ylim([-10 10])


x_total=x_rec2;

figure(13)
%plot(t(1:row_x2),x_total,'k',t(1:row_x2),x(1:row_x2),'r:')
plot(t,x_total,'k',t,x,'r:')
xlabel('Time/s');
ylabel('x');
legend('Reconstructed','Original');

figure(14)
%plot(t(1:row_x2),x_remain(:,1),'b')
plot(t,x_remain,'k',t,z2,'r:')
xlabel('Time/s');
ylabel('x_2');
legend('Reconstructed','Original');

%-------------------------------------------------------
figure(15)
%subplot(5,1,1),plot(t,x_total,'k',t,x,'r:')
%subplot(5,1,2),plot(t,x_rec2(:,1),'k',t,z1,'r:')
%subplot(5,1,3),plot(t,x_rec2(:,2),'k',t,z2,'r:')
%subplot(5,1,4),plot(t,x_rec2(:,3),'k',t,z3,'r:')
%subplot(5,1,5),plot(t,x_remain,'k')

subplot(3,1,1);plot(t,x,'k');ylabel('x');set(gca,'XTickLabel',[])
subplot(3,1,2);plot(t,x_rec2(:,1),'k');ylabel('H_1');set(gca,'XTickLabel',[])
subplot(3,1,3);plot(t,x_remain,'k');ylabel('H_r');
xlabel('Time/s');


f=0:1/(num_t*Ts):1/Ts;
x_f1=fft(x_rec2(:,1));
x_f2=fft(x_remain);

figure(16)
plot(f(1:(numel(f)+1)/2),abs(x_f1(1:(numel(f)+1)/2)),'k');ylabel('H_1')
%subplot(2,1,1);plot(f(1:(numel(f)+1)/2),abs(x_f1(1:(numel(f)+1)/2)),'k');ylabel('H_1');set(gca,'XTickLabel',[])
%subplot(2,1,2);plot(f(1:(numel(f)+1)/2),abs(x_f2(1:(numel(f)+1)/2)),'k');ylabel('H_2');

xlabel('Frequency/Hz');