SAM_SSD.m

function [varargout]= SAM_SSD（varargin）
%例 SSC=SAM_SSD（yfs）;
%分解结果不包含残余分量
%输入：信号，采样频率，最大分解次数
warning off; %关闭警告

if nargin==2
    y=varargin{1};
    Fs=varargin{2};
    Times=inf;
    Nstdd=1.2;
elseif nargin==3
    y=varargin{1};
    Fs=varargin{2};
    Times=varargin{3};
    Nstdd=1.2;
elseif nargin==4
    y=varargin{1}; %原始信号
    Fs=varargin{2}; %采样频率
    Times=varargin{3}; %输出SSC的数量（默认不设置）
    Nstdd=varargin{4}; %嵌入维数权值（默认0.2）
end

%输入： 待分解的时域信号y 信号采样频率  
    k1=0;
    th=0.01; %迭代停止阈值
    y=y（:）‘; %强制变为行向量
    L=length（y）;

    y = y-mean（y）; %将输入信号变为均值为0的向量，因为SSD针对的就是均值为0的向量

    if Fs/L <= 0.5
        lf = L;
    else
        lf = 2*Fs;
    end
    RR1=zeros（size（y））; %初始化输出
    orig = y;
    remen = 1;  %迭代停止标准
    testcond = 0;
    time=0;
    
    while （remen > th&&time        time=time+1;
        k1 = k1+1;
        y = y-mean（y）;  %输入信号均值向0逼近
        v2=y;
        clear nmmt;
        [nmmtff] = pwelch（v2[][]4096Fs）;%求功率谱估计密度
        [~in3] = max（nmmt）;
        nmmt = nmmt‘;
        if （（k1 == 1） && （ff（in3）/Fs < 1e-3）） % 在第一次迭代中，如果最大PDS值/采样频率<阈值（10^-3），则认为残余量具有单调趋势，将嵌入维度M设置为N/3  否则 M=1.2*FS/最大PDS值
            l = floor（L/3）;     %l就是嵌入维数    %设定嵌入维度，可更改的地方
            %下面开始重构SSC分量
            %先是重构轨迹矩阵
            M = zeros（L-（l-1） l）;
            for k=1:L-（l-1）
                M（k:） = v2（k:k+（l-1））;
            end
            [USV] = svd（M）;   %对M进行奇异值分解 M=U*S*V   U左奇异向量 S奇异值 V右奇异向量
            U（:l+1:end） = [];
            S（l+1:end:） = [];
            V（:l+1:end） = [];

            rM = rot90（U（:1）*S（1:）*V‘）;  %将数组旋转90°
            r = zeros（1L）;
            [~m] = size（rM）;
            for k=-（l-1）:L-（l）
                r（k+l） = sum（diag（rMk））/m;  
            end
        else % 以后的运行中

            for cont = 1:2 %运行2次
                v2 = v2-mean（v2）;
                [deltaf] = gaussfit_SSD（ffnmmt‘）; % 高斯拟合
                % 带宽估计
                [~iiii1] = min（abs（ff-（ff（in3）-deltaf）））;
                [~iiii2] = min（abs（ff-（ff（in3）+deltaf）））;
                
                 l = floor（Fs/ff（in3）*Nstdd）;   %设定嵌入维度，可更改的地方
       
                if l <= 2 
                    l=2;
                elseif l > floor（L/3）
                    l = floor（L/3）;
                end
                
                M=zeros（L l）;
                % M built with wrap-around
                for k=1:l
                    M（:k）=[v2（end-k+2:end）‘; v2（1:end-k+1）‘];
                end
                [USV] = svd（M0）;

                %选择主成分
                 if size（U2）>l
                    yy = abs（fft（U（:1:l）lf））;
                else
                    yy = abs（fft（Ulf））;
                end
                yy_n = size（yy1）;
                ff2 = （0:yy_n-1）*Fs/yy_n;
                yy（floor（yy_n/2）+1:end:） = [];
                ff2（floor（length（ff2）/2）+1:end） = [];
                % %