奇異譜分析

步驟一：建立軌跡矩陣

原始信號長度為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