奇异谱分析

步骤一：建立轨迹矩阵

原始信号长度为N，滑动窗口长度为Lp，Kp = N-Lp+1；轨迹矩阵就是按照列做分割，第一列为索引为1~Lp的信号，第二列为2~Lp+1，第三列为3~Lp+2，第Kp列为信号索引为Kp~N。

 轨迹矩阵：

 

步骤二：奇异值分解

1） 计算XX^T的特征值和特征向量U

2）  计算左奇异向量U和右奇异向量V，

求V的时候可以不用除lambda，因为重构信号的时候又乘上lambda。

步骤三：分组

分组的目的就是将目标信号成份和其他信号成份分开，在信号处理领域，通常认为前面r个较大的奇异值反应信号的主要能量。

步骤四：对角重构信号平均化

根据分组结果将对应的奇异向量重构：

i为选择的r个奇异向量。

对角平均化分为三部完成，对应于下面表格的三部分。

若：奇异矩阵是rca，Lp*Kp，其中Lp<Kp，重构信号为y，长度为N

第一部分：浅蓝色部分，1~Lp-1

y(1) = rca(1,1);

y(2) = (rca(1,2)+rca(2,1))/2;

y(3) = (rca(1,3)+rca(2,2)+rca(3,1))/3;

…

y(Lp-1) = (rca(1,Lp-1)+rca(2,Lp-2)+…+rca(Lp-1,1))/(Lp-1);

第二部分：橙色部分，Lp~Kp

y(Lp) = (rca(1,Lp)+rca(2,Lp-1)+…+rca(Lp,1))/Lp;

y(Lp+1) = (rca(1,Lp+1)+rca(2,Lp)+rca(3,Lp-1)…+rca(Lp,2))/Lp;

…

y(Kp) = (rca(1,Kp)+rca(2,Kp-1)+rca(3,Kp-2)+…+rca(Lp,Kp-Lp+1))/Lp;

第三部分：绿色部分，Kp+1~N

y(Kp+1) = (rca(2,Kp)+rca(3,Kp-1)+rca(4,Kp-2)+…+rca(Lp, Kp-Lp+2))/(Lp-1);

y(Kp+2) = (rca(3,Kp)+rca(4,Kp-1)+…+rca(Lp,Kp-Lp+3))/(Lp-2)

…

y(N-1) = (rca(Lp-1,Kp)+rca(Lp,Kp-1))/(Lp-(Lp-1)+1);

y(N) = rca(Lp,Kp);

function [signalFiltered]=SSA(signal,windowLen)
%========================================================================
% 参看http://www.ilovematlab.cn/thread-301868-1-1.html柳絮艳的分享代码ssa改写
% signal 原始信号
% windowLen 窗口长度
% signalFiltered   重构时间序列
%=========================================================================
% Step1 : 建立轨迹矩阵
N=length(signal);
if windowLen>N/2;
    windowLen=N-windowLen;
end
K=N-windowLen+1;
X=zeros(windowLen,K);
for i=1:K
    X(1:windowLen,i)=signal(i:windowLen+i-1);
end
% Step 2: 奇异值分解
S=X*X';
[U,autoval]=eig(S);%eig返回矩阵的特征值和特征向量，U是特征向量，autoval是特征值
[d,i]=sort(diag(autoval),'descend');
sev=sum(d); %特征值的和,可根据占比选择有效信号
U=U(:,i);
V=(X')*U;
% Step 3:分组
I=[1:windowLen/2];%I的选择可根据信号特征选择
Vt=V';
rca=U(:,I)*Vt(I,:);
% Step 4: 对交平均化重构信号
y=zeros(N,1);
Lp=min(windowLen,K);
Kp=max(windowLen,K);
%重构 1~Lp-1
for k=0:Lp-2
    for m=1:k+1;
        y(k+1)=y(k+1)+(1/(k+1))*rca(m,k-m+2);
    end
end
%重构 Lp~Kp
for k=Lp-1:Kp-1
    for m=1:Lp;
        y(k+1)=y(k+1)+(1/(Lp))*rca(m,k-m+2);
    end
end
%重构 Kp+1~N
for k=Kp:N-1
    for m=k-Kp+2:N-Kp+1;
        y(k+1)=y(k+1)+(1/(N-k))*rca(m,k-m+2);
    end
end

signalFiltered = y;
end

　　

　　

参考：[1] https://wiki.mbalib.com/wiki/%E5%A5%87%E5%BC%82%E8%B0%B1%E5%88%86%E6%9E%90

[2]《基于改进奇异谱分析的信号去噪方法》http://journal.bit.edu.cn/zr/ch/reader/create_pdf.aspx?file_no=20160713&year_id=2016&quarter_id=7&falg=1