function [Upi,xf,V] = ssa(ts,L),
N=length(ts);
K=N-L+1; %K受制于ts的长度N和滑动窗口L
X=zeros(L,K);
for i=1:L, %构造轨迹矩阵L行K列，原理是从第一个开始按ts顺序取K个当第一行，再从第二个开始取K个当第二行，依次取到L行
X(i,:)=ts(i:i-1+K)';
end;
S=X*transpose(X);
[U,Lemda]=eig(S); %eig(S)指对s进行对角化，含特征向量的矩阵U，Lemda是对角矩阵
lemda=diag(Lemda); %把lemda元素排成列向量形式 这里的diag（相量）等于对角，diag（对角）=相量
[Y,I]=sort(lemda,'descend'); %对列向量lemda做降序排列赋给Y，并且I作为其索引，I(i)表示降序后的第i个元素在原lamda的位置，见同文件的另一个示例
for k=1:length(Y),
if lemda(I(k))>0, %表示从大到小检索Y的元素（S的特征值），当等于0时停止循环，aa'矩阵的性质是特征值必大于等于0，并且可以把0滤了
Upi(:,k)=U(:,I(k));%使Upi矩阵的列向量按特征值从大到小对应的特征向量排列，如Upi(:,1)=U(:,I(1))即将特征值按降序排列时的第一个元素（也就是特征值最大那个逼）在原位置的序号，对应的U的那个列赋给Upi的第一列
end;
end;
D=size(Upi,2) %函数size(A,1)、size(A,2)和size(A,3)返回矩阵对应的行数、列数、维度的值
%此时做了Upi成为单位正交矩阵，若eig对象为对称，则自动生成单位正交的特征向量拼成的矩阵
for k=1:D,% D意思等价于特征值不为零的个数
V(:,k)=transpose(X)*Upi(:,k)/sqrt(lemda(I(k))); %算出X'X对应的单位正交的特征向量并按列拼成矩阵V
Xpi(k,:,:)=sqrt(lemda(I(k)))*Upi(:,k)*transpose(V(:,k));% Xpi(k,:,:)相当于是Xpi包含三维，第k个矩阵的：行，：列，整句的意思就是把X矩阵化成L个（奇异值*对应那两个特征向量）赋给Xpi
end;
for d=1:D,
Xtutor=zeros(L,K);
for k=1:K,
Xtutor(:,k)=Xpi(d,:,k);%意思是Xpi中第d个矩阵的第k列。一整局的意思就是这个for完成后的Xtutor就是（奇异值*对应那两个特征向量）
end;
xf(d,:)=Hankelization(Xtutor,L,K,N);%这个矩阵相当于是把Xtutor对角平均，并且赋给xf的每一行
end;
end
评论0