计算原始信号减去EEMD（Ensemble Empirical Mode Decomposition，集成经验模态分解）、SSA（Singular Spectrum Analysis，奇异谱分析）和LS（Least Squares，最小二乘法）季节模型所产生的残差的频谱指数的matlab代码

以下是计算原始信号减去EEMD、SSA和LS季节模型所产生的残差的频谱指数的matlab代码：

% 假设原始信号为signal，采样率为fs
% EEMD分解
imf = ceemdan(signal); % 进行CEEMDAN分解，得到各个本征模态函数imf
residual = signal - sum(imf, 2); % 计算残差
% 计算频谱指数
[pxx, f] = pwelch(residual, [], [], [], fs); % 计算残差的功率谱密度pxx和对应频率f
fi = find(f > 0 & f < 0.1); % 选择感兴趣的频率区间
SI_EEMD = trapz(pxx(fi)) / trapz(pxx); % 计算频谱指数SI_EEMD

% SSA分解
L = 100; % 选择时间窗口长度
K = floor(length(signal) / L); % 计算分解后的子序列数
X = zeros(L, K);
for i = 1:K
    X(:, i) = signal((i-1)*L+1:i*L);
end
imf = ssa(X); % 进行SSA分解，得到各个本征模态函数imf
residual = signal - sum(imf, 2); % 计算残差
% 计算频谱指数
[pxx, f] = pwelch(residual, [], [], [], fs); % 计算残差的功率谱密度pxx和对应频率f
fi = find(f > 0 & f < 0.1); % 选择感兴趣的频率区间
SI_SSA = trapz(pxx(fi)) / trapz(pxx); % 计算频谱指数SI_SSA

% LS季节模型
t = (1:length(signal))'; % 构造时间序列t
X = [ones(length(signal), 1), sin(2*pi*t/365), cos(2*pi*t/365)]; % 构造基函数矩阵X
beta = X \ signal; % 进行最小二乘拟合，得到拟合系数beta
residual = signal - X*beta; % 计算残差
% 计算频谱指数
[pxx, f] = pwelch(residual, [], [], [], fs); % 计算残差的功率谱密度pxx和对应频率f
fi = find(f > 0 & f < 0.1); % 选择感兴趣的频率区间
SI_LS = trapz(pxx(fi)) / trapz(pxx); % 计算频谱指数SI_LS

其中，ceemdan和ssa是EEMD和SSA的matlab函数，可在MATLAB File Exchange中下载。