奇异谱分析

步骤一：建立轨迹矩阵

原始信号长度为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

y1) = rca1,1);

y2) = rca1,2)+rca2,1))/2;

y3) = rca1,3)+rca2,2)+rca3,1))/3;

…

yLp-1) = rca1,Lp-1)+rca2,Lp-2)+…+rcaLp-1,1))/Lp-1);

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

yLp) = rca1,Lp)+rca2,Lp-1)+…+rcaLp,1))/Lp;

yLp+1) = rca1,Lp+1)+rca2,Lp)+rca3,Lp-1)…+rcaLp,2))/Lp;

…

yKp) = rca1,Kp)+rca2,Kp-1)+rca3,Kp-2)+…+rcaLp,Kp-Lp+1))/Lp;

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

yKp+1) = rca2,Kp)+rca3,Kp-1)+rca4,Kp-2)+…+rcaLp, Kp-Lp+2))/Lp-1);

yKp+2) = rca3,Kp)+rca4,Kp-1)+…+rcaLp,Kp-Lp+3))/Lp-2)

…

yN-1) = rcaLp-1,Kp)+rcaLp,Kp-1))/Lp-Lp-1)+1);

yN) = rcaLp,Kp);

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