奇异谱分析（SSA）matlab

最近在做的是通过将一个正弦函数（三个不同频率），分解重构，原理过程可以参考下面两个

[机器学习] 奇异谱分析(SSA)原理及Python实现-CSDN博客

奇异谱分析(SSA)的matlab实现_奇异谱分析 matlab-CSDN博客

运行的文件是SSA1.m，其他的.m文件都是SSA调用的。具体代码如下：

SSA1.m

%输入合成信号y
f1=50;
f2=90;
f3=410;
fs=50;%采样频率
t=0:1/fs:1;
y=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);
plot(t,y,'r*-');
%ssa
n = length(y); %测试数据长度
l = 10;        %输出窗口长度L
[y1,lamuda]=SSA_function(y,n,l);         %进行奇异谱分析得到L个初等矩阵
ref = 0.95;
[x] = Contribution_rate(lamuda,l,ref);   %根据贡献率选择重构阶数，阈值根据需求设置
[coll] = w_collelation(y1,l,n);          %权相关分析得到w-correlation图
figure(1)
HeatMap(coll);
yout = zeros(1,n);
for i=1:l
    subplot(5,2,i);
    plot(y1(i,:))
end    
for i=1:x
   yout(1,:) = yout(1,:) + y1(i,:); 
end
figure(2)
plot(1:n,y,'-r',1:n,yout,'-g')
%% 输出趋势项
prompt = '趋势项阶数: ';              %根据w-collelation图输入重构阶数
x = input(prompt);                       %根据权重图选择要重构信号
yout = zeros(1,n);
for i=1:x
   yout(1,:) = yout(1,:) + y1(i,:); 
end
figure(3)
plot(1:n,yout,'-r')
%% 输出周期项
prompt = '输出周期项: ';              %根据w-collelation图输入重构阶数
x2 = input(prompt);                       %根据权重图选择要重构信号
yout = zeros(1,n);
for i=x+1:x2
   yout(1,:) = yout(1,:) + y1(i,:); 
end
figure(4)
plot(1:n,yout,'-g')

SSA_function.m

%% 奇异谱分析函数
%x 原始时间序列；n一维时间序列x的维度；l窗口长度(为周期的整数倍，不超过n/2)
% 返回的y的维度是[窗口长度l*序列长度n],表示l个初等矩阵；
function [y,lamuda]=SSA_function(x,n,l)
    % Step1 : 建立时滞矩阵
    k1=n-l+1;
    X2=zeros(l,k1);
    for i=1:k1
        X2(1:l,i)=x(i:l+i-1);
    end
    % Step 2: 奇异值分解
    x3=X2*X2';
    [U,S,V] = svd(x3);
    for i=1:l
        V1(:,i) = X2'*U(:,i)/sqrt(S(i,i)); 
    end
    %初等矩阵重构得到l个时间序列；
    y = zeros(l,n);
    Lp=min(l,k1);
    Kp=max(l,k1);
    for i=1:l
        xi = sqrt(S(i,i)) * U(:,i) * V1(:,i)';
        y1=zeros(n,1);
        for k=0:Lp-2
            for m=1:k+1
                y1(k+1)=y1(k+1)+(1/(k+1))*xi(m,k-m+2);
            end
        end
        %重构 Lp~Kp
        for k=Lp-1:Kp-1
            for m=1:Lp
                y1(k+1)=y1(k+1)+(1/(Lp))*xi(m,k-m+2);
            end
        end
        %重构 Kp+1~N
        for k=Kp:n-1
            for m=k-Kp+2:n-Kp+1
                y1(k+1)=y1(k+1)+(1/(n-k))*xi(m,k-m+2);
            end
        end 
        y(i,:) = y1(:,1);
        
            
    end
    %将特征值输出一列
    %lamuda = zeros(l,1)
    for i = 1:l 
        lamuda(i,1) = S(i,i);
    end
end

w_collection.m

%% 权相关
function [coll] = w_collelation(y,l,n)
    k = n-l+1;
    coll = zeros(l,l);
    for i =1:l
        for j =1:l
            if l<k
                xfc = 0.0;
                xfx = 0.d0;
                xfy = 0.0;
                for h=1:n
                    if h >= 1 && h<l
                        w = h;
                    elseif h >= l && h<k
                        w=l;
                    elseif h >= k && h<=n
                        w = n-h+1;
                    end
                    xfc = xfc + w*y(i,h)*y(j,h);
                    xfx = xfx + w *y(i,h)^2;
                    xfy = xfy + w * y(j,h)^2;
                end
                coll(i,j) = xfc/(sqrt(xfx)*sqrt(xfy));
            end
        end
    end
end     

Contrabution_rate.m

%% 计算奇异值的贡献率
function [n,array2] = Contribution_rate(lamuda,l,ref) %输入奇异值、窗口长度、阈值(0,1)
    sum1 = sum(lamuda(:));
    for i=1:l
        sum2 = sum(lamuda(1:i));
        num = sum2/sum1;
        array2(i) = num;
        if num >= ref
            n = i;
            break
        end
    end
end   

实现效果如下：

heatmap是原始图与重构图的相关性，可以根据权重来改变参数。最右边是不同奇异值对应的子图，左下图是原始与重构的图，左上是原始图。