奇异谱分析 SSA matlab

参考链接：https://www.kaggle.com/code/jdarcy/introducing-ssa-for-time-series-decomposition/notebook

%% 数据
clear
clc
close  all
t = [0:200];
p1 =20;
p2 =30;
f1 = 0.001*(t-100).^2 ;
f2 = 2*sin(2*pi*t./p1);
f3 = 0.75*sin(2*pi*t./p2) ;
f4 = rand(size(t))*2 -1;
f = f1+f2 +f3+f4;
figure;
plot(f)
hold on
plot(f1);
hold on
plot(f2 );
hold on
plot(f3);
hold on
plot(f4);
xlim([0 numel(t)]);
legend('结果：f','抛物线：f1','正弦1：f2','正弦2：f3','噪声：f4')
%% 轨迹矩阵 X
L =70;%L为窗口长度，窗口内的数据为一行，然后以步长为1，将窗口向右移动，直至窗口覆盖f的最后一个数据

if L<2 | L>numel(f)/2
    error('窗口长度不合理');
end
N =numel(f);%数据的长度
K = N - L+1; % 窗口可滑动次数, X 的列
X = zeros(L,K);% L 为行数，窗口长度的值为行数，K为列数，每一列的数据为一个窗口内的数据 ，
% X1表示 X矩阵 第一列的数据 (ft1...ftL)‘  其中’表示转置 。其中共有 K 列

for i =1:K % 获取 X矩阵 每一列应该的数据  
    Xn = f(i:i+L-1)';  % X矩阵中每一列的数据为 窗口从1开始向后移动中窗口内的数据 
    X(:,i) = Xn; % 第k次移动窗口后，窗口内的数据，赋给X矩阵的第k列
end

%% 对轨迹矩阵X 进行奇异值分解
[u,a,v] = svd(X); 
% 直接调用了奇异值分解函数 
% u 为 L*L 的酉矩阵，包含了X的左奇异向量正交集合作为列
% a 为 L*K 的对角矩阵， 包含了X的 n 个奇异值，从大到小排列；
% k 为 K*K 的酉矩阵， 包含了X的右奇异向量正交集合作为列

% 轨迹矩阵的SVD分解结果重写
% 其中重写 X 矩阵的过程为： 
% a矩阵中选择一个奇异值，
% 用该奇异值对应的行号，从 u 矩阵中选出对应列号的一列 (L*1)
% 用该奇异值对应的列号，从 v 矩阵中选出对应列号一列 (K*1)的转置v'(1*K)
% 将 u 矩阵中选出对应列号的一列 (L*1)， 与v 矩阵中选出对应列号一列 (K*1)的转置v'(1*K)相乘
% 然后将该矩阵（L*K）与对应的奇异值相乘，得到对应的矩阵a*（L*K）
% 然后将上述步骤，每一个奇异值所得矩阵a*（L*K）累加，得到的即为重写后 X
% 把前 r < n个  上述每列所得 初等矩阵a*（L*K）加起来，即可得一个 X 矩阵的低维近似。

% close all

X1 = zeros(L,K); % 用于存储新 X 矩阵 
Yi = zeros(L,L+K-1); % 用于存储每个初等矩阵转换为一维度数据后的数据
X2 = zeros(L,K,L+K-1);% 用于存储每个 对线线平均化后的 初等矩阵 一共L+K-1个 
X3 = zeros(L,K,L+K-1);% 用于存储每个 对线线平均化后的 初等矩阵 一共L+K-1个 
X5 = zeros(L,K,L+K-1);% 用于存储每个 对线线平均化后的 初等矩阵 一共L+K-1个

for i=1:L %全部显示的话就是 L 遍历矩阵中的每一行的数据 ，由于带有显示所以用了 11代替L
% for i=7:10
%     for j =1:k % 注释掉K是因为 a矩阵中 L:end 全为零 ，没有计算意义，所以博客写的都是 ai*ui*v'i
    Xi = a(i,i)*u(:,i)*v(:,i)';
    X1 = X1+Xi;
%     figure;
%     imshow(Xi,[]) % 用于显示初等矩阵
%     title([' 第 ' num2str(i) '个初等矩阵示意图']);
%     colormap('hot')

    %对角线平均 还原成一维数据  直接写成一个函数吧 全放这就太乱了 
%     figure;
    [Yi(i,:),X2(:,:,i),X3(:,:,i),X5(:,:,i)] = djx1(Xi);
%     plot(Yi(i,:));  % 用于显示初等矩阵对应的一维数据
%     
%     figure; % 这个是对角线单个数据的 是没有意义 写错的 
%     imshow( X2(:,:,i) ,[]) % 用于初等矩阵 对角线 单 个 数据除以对角线个数值后的 矩阵 
% %     imshow(reshape(mapminmax(reshape(a,1,15)),5,3) ,[]) %  用做了归一化 但是没啥意义 因为显示的时候加了 []
%     title([' 第 ' num2str(i) '个初等矩阵 单个 变换示意图']); 
%     colormap('hot') 
    
%     figure; % 这个是对角线累加数据的 可以用 跟参考的不一样 但是能用
%     imshow( X3(:,:,i) ,[]) % 用于初等矩阵 对角线数据 累 加 值 除以对角线个数值后的 矩阵 
% %     imshow(reshape(mapminmax(reshape(a,1,15)),5,3) ,[]) %  用做了归一化 但是没啥意义 因为显示的时候加了 []
%     title([' 第 ' num2str(i) '个初等矩阵 累加 变换示意图']);
%     colormap('hot')
%     
%     figure; % 这个是对角线累加数据的 可以用 跟参考的不一样 但是能用
%     imshow( X5(:,:,i) ,[]) % 用于初等矩阵 对角线数据 累 加 值 除以对角线个数值后的 矩阵 
% %     imshow(reshape(mapminmax(reshape(a,1,15)),5,3) ,[]) %  用做了归一化 但是没啥意义 因为显示的时候加了 []
%     title([' 第 ' num2str(i) '个初等矩阵 累加 变换示意图']);
%     colormap('hot')

end
sum(sum(abs(X1- X))) %累加 上述 i 个初等矩阵 时后，与原图的重写误差
%%  单个Xi矩阵重构  对角线平均化 
% 观察曲线后 叠加一起   对角线数据除以个数值后的初等矩阵没观察出来
figure;
plot(Yi(1,:)+Yi(2,:)+Yi(7,:)); 
hold on
plot(Yi(3,:)+Yi(4,:)); % 3和4 与 5和6的周期不一样
hold on
plot(Yi(5,:)+Yi(6,:));
hold on
plot(sum(Yi(8:11,:)));

效果图：

对角线平均函数(单独建立一个djx1.m)：

function [temp,Xi,X4,X5] = djx1(Xi)
[L,K] = size(Xi);
m=0;
if L > K
    temp = K;
    K = L;
    L =temp;
    Xi =Xi';
    m =1;
end
fi =zeros(1,L+K-1);% 累积对角线的值
fit = zeros(1,L+K-1); %对对角线值个数计数
X3 = zeros(L,K);%对角线累加矩阵 
X4 = zeros(L,K);%对角线累加矩阵 
X3(1,:) = Xi(1,:);
X3(:,K) =  Xi(:,K);
X4(1,:) = Xi(1,:);
X4(:,K) =  Xi(:,K);
for i=1:L  %遍历Xi矩阵所有行数据
    for j=1:K % 遍历Xi矩阵所有列数据
            fi(1,i+j-1) = fi(1,i+j-1)+Xi(i,j);
            fit(1,i+j-1) = fit(1,i+j-1)+1;
            if  fit(1,i+j-1)>1 & j+1<=K
                X3(i,j) =Xi(i,j)+ X3(i-1,j+1) ; % 当前X3(i,j) 为空  X3(i-1,j-1)为上一个累加值  Xi(i,j) 为此处应该的值
                X4(i,j) = X3(i,j)/ fit(1,i+j-1) ;
            end
    end
end
temp = fi./fit;
% fi
% fit
% temp

for i=1:L  %遍历Xi矩阵所有行数据
    for j=1:K % 遍历Xi矩阵所有列数据
            Xi(i,j) =Xi(i,j)/fit(1,i+j-1) ;
    end
end

% 上个步骤只是矩阵中 单个数据的均值 而矩阵转化为一位序列时 为累加的均值 所以如果想通过图观测曲线则应该用累积的均值

% 把曲线数据所对应的均值填回去
X5 = zeros(L,K+1);
for i=1:L  %遍历Xi矩阵所有行数据
    for j=1:K % 遍历Xi矩阵所有列数据
        X5(i,j) = temp(1,i+j-1);
    end
end
X5(:,K+1) = []; 
if m == 1
Xi = Xi';
X4 = X4';
X5 = X5';
end

end