奇异谱分析

Posted 2021-01-26 xhslovecx

步骤一：建立轨迹矩阵

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

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

[2]《基于改进奇异谱分析的信号去噪方法》