ssa_1.m

一个关于混沌序列的噪声压缩程序

%Code for decomposing univariate time series via
%Caterpillar-SSA method
%
%The programme first folds the univariate series
%into mulitvariate ones, thus creating m new
%time series. Next, the covariance matrix of these
%series is calculated, which is then decomposed via
%SVD. The folded series is projected onto the space
%spanned by its principal components and reconstructed
%in terms of the original time series.
clear all
close all
% load('c:\MATLAB7\work\noisechain\x1.mat');      % 装入滤波后的时间序列
% x1=x0;
% x1=x0(:,1)

disp('---------- Lorenz ----------')

sigma = 16;             % Lorenz 方程参数 a
b = 4;                  %                 b
r = 45.92;              %                 c             

y = [-1,0,1];           % 起始点 (1 x 3 的行向量)
h = 0.01;               % 积分时间步长

k1 = 1000;             % 前面的迭代点数
k2 = 200;              % 后面的迭代点数


k1 = 10000;             % 前面的迭代点数
k2 = 500;               % 后面的迭代点数

z = LorenzData(y,h,k1+k2,sigma,r,b);
x = z(k1+1:end,1);      % 时间序列
figure
plot(x)


% Loranz+noise 加200％的噪声？？？？
N=length(x);
wn=20*randn(N,1);
x1=x+wn;

figure
plot(x1)
%Creating (k by m) matrix from the original series
%m=length of Caterpillar, approximation to the true
%dimension of the process

	T=length(x1);
	m=7;
	k=T-m+1;
	for i=1:m
	xx1(1:k,i)=x1(i:k+i-1);
	end

%Normalizing the observations matrix:
%Subtracting column means and dividing by column
%standard deviations

for i=1:m
   xavr(i)=mean(xx1(:,i));
   xdev(i)=std(xx1(:,i));
end
xavr=ones(k,1)*xavr;
xx1=xx1-xavr;

for i=1:m
   for j=1:k
    xx1(j,i)=xx1(j,i)/xdev(i);
   end
end

%Eigenvalue-eigenvector decomposition of covariance matrix
% 	z=xx1'*xx1;
% 	[v,d]=eig(z);
% 	dd=diag(d);
%     n=length(dd);
    
kernel='gaussian';
kerneloption=1;
%max_eigvec=8;
[v,d]=kernelpca(xx1,kernel,kerneloption);
dd=diag(d);
n=length(dd);
    
    
 for i=1:n-1                         %特征值从大到小排序，特征向量也相应调整次序
    for j=i+1:n
        if dd(i)<dd(j)
            temp=dd(i);
            temp1=v(:,i);
            dd(i)=dd(j);
            v(:,i)=v(:,j);
            dd(j)=temp;
            v(:,j)=temp1;
        end
    end
 end
dd=dd./sum(dd);
l=find(dd>0.3);   %提取特征值大于1的特征向量
v=v(:,l);

figure
plot(dd)
title('Eigenvalues of Autocovariance Matrix')


%Computing principal components and
%reconstructing original time series
%One can restore any subset of PCs: from 1 to m
pc=xx1*v;
rc=pc*v';

%Denormalizing the RCs
rc=rc+xavr;
for i=1:m
   for j=1:k
    rc(j,i)=rc(j,i)*xdev(i);
   end
end

%重构时间序列

y=zeros(T,1);
for i=1:k
  for j=1:m
     if i<=k
	 y(i+j-1,1)=y(i+j-1,1)+rc(i,j);
     else
	 y(i+j,1)=y(i+i)+rc(i,j);
     end
   end
end

for i=1:T
     if i<=m
	  y(i,1)=(y(i,1)/i);
     elseif i<=T-m
	  y(i,1)=(y(i,1)/m);
     else
	 y(i,1)=(y(i,1)/(T-i+1));
     end
end
y=y-mean(y);
figure
plot(y)
r=sum((x-y).^2)/sum((x1-x).^2);

%计算奇异谱

s1=svd(xx1);
s1=log(s1./sum(s1));
s2=svd(rc);
s2=log(s2./sum(s2));
figure
plot(s1)
figure
plot(s2)