function [num,Period, Frequency, Density, CL95]=spectrum(x,mLAG) %%% function for power spectral analysis % usage: [num,Period, Frequency, density, cl95]=spectrum(x,mLAG) % Gong Daoyi 2003.12 xLEN=length(x); SER=x;N=xLEN;mLAGWK=mLAG;mLEN=N;J=mLAG;J1=J+1; %c calculating auto-connection coefficient A=0.0; C=0.; for I=1:N A=A+SER(I);end % I A=A/N; for I=1:N SER(I)=SER(I)-A; C=C+SER(I).^2; end % I C=C/N; for L=1:J CC(L)=0.0; for I=1:N-L CC(L)=CC(L)+SER(I)*SER(I+L); end %I CC(L)=CC(L)./(N-L); CC(L)=CC(L)/C; end %L C=1.0; %c estimating rude power spectra SPE(1)=0.0; for L=1:J-1 SPE(1)=SPE(1)+CC(L); end %L SPE(1)=SPE(1)./J+(C+CC(J))./(2*J); for L=1:J-1 % DO 210 L=1,J-1 SPE(L+1)=0.; for I=1:J-1 SPE(L+1)=SPE(L+1)+CC(I)*cos(pi*L*I/J); end % I SPE(L+1)=2*SPE(L+1)./J+C./J+(-1).^L*CC(J)./J; end % 210 L SPE(J1)=0.0; for I=1:J-1 SPE(J1)=SPE(J1)+(-1).^I*CC(J); end %I SPE(J1)=SPE(J1)/J+(C+(-1).^J*CC(J))/(2*J); %c smoothing power spectra PS(1)=.54*SPE(1)+.46*SPE(2); for L=2:J PS(L)=.23*SPE(L-1)+.54*SPE(L)+.23*SPE(L+1); end %L PS(J1)=.46*SPE(J)+.54*SPE(J1); %c statistical significence of PS W=0.0; for L=1:J-1 W=W+SPE(L+1); end %L W=W/J+(SPE(1)+SPE(J1))/(2*J); if (J > fix(N/2)) W=2.57*W; end if(J == fix(N/2)) W=2.49*W; end if(J < fix(N/2) & J > fix(N/3)) W=2.323*W; end if (J == fix(N/3)) W=2.157*W; end if (J < fix(N/3)) W=1.979*W; end %c the red noice examination for L=1:J1 SK(L)=W*(1-CC(1).^2)/(1+CC(1).^2-2*CC(1)*cos(3.14159*(L-1)/J)); end % L if (CC(1) > 0 & CC(1) >= CC(2) ) %c the white noice examination else for L=1:J1 SK(L)=W; end %L end % if %c calculating the length of cycle T(1)=NaN; for L=2:J1 T(L)=(2.0*J)/(L*1.0-1.0); end % L num=1:J+1;num=num(:)-1; Period=T(:); Frequency=1./T(:); Density=PS(:); CL95=SK(:);

这段代码是MATLAB中进行功率谱分析的函数spectrum，用于分析时间序列信号的频谱特性。下面是对代码的详细解释：

function [num,Period, Frequency, Density, CL95]=spectrum(x,mLAG)
% x: 时间序列信号
% mLAG: 滞后阶数

xLEN=length(x); % 获取信号长度
SER=x; % 将信号存储到SER数组中
N=xLEN;
mLAGWK=mLAG;
mLEN=N;
J=mLAG;
J1=J+1;

% 计算自相关系数
A=0.0;
C=0.;
for I=1:N
    A=A+SER(I);
end
A=A/N;
for I=1:N
    SER(I)=SER(I)-A;
    C=C+SER(I).^2;
end
C=C/N;
for L=1:J
    CC(L)=0.0;
    for I=1:N-L
        CC(L)=CC(L)+SER(I)*SER(I+L);
    end
    CC(L)=CC(L)./(N-L);
    CC(L)=CC(L)/C;
end

C=1.0;

% 估计粗略的功率谱
SPE(1)=0.0;
for L=1:J-1
    SPE(1)=SPE(1)+CC(L);
end
SPE(1)=SPE(1)./J+(C+CC(J))./(2*J);
for L=1:J-1
    SPE(L+1)=0.;
    for I=1:J-1
        SPE(L+1)=SPE(L+1)+CC(I)*cos(pi*L*I/J);
    end
    SPE(L+1)=2*SPE(L+1)./J+C./J+(-1).^L*CC(J)./J;
end
SPE(J1)=0.0;
for I=1:J-1
    SPE(J1)=SPE(J1)+(-1).^I*CC(J);
end
SPE(J1)=SPE(J1)/J+(C+(-1).^J*CC(J))/(2*J);

% 平滑功率谱
PS(1)=.54*SPE(1)+.46*SPE(2);
for L=2:J
    PS(L)=.23*SPE(L-1)+.54*SPE(L)+.23*SPE(L+1);
end
PS(J1)=.46*SPE(J)+.54*SPE(J1);

% 计算功率谱的95%置信区间
W=0.0;
for L=1:J-1
    W=W+SPE(L+1);
end
W=W/J+(SPE(1)+SPE(J1))/(2*J);
if (J > fix(N/2))
    W=2.57*W;
end
if(J == fix(N/2))
    W=2.49*W;
end
if(J < fix(N/2) & J > fix(N/3))
    W=2.323*W;
end
if (J == fix(N/3))
    W=2.157*W;
end
if (J < fix(N/3))
    W=1.979*W;
end

% 进行红噪声检验
for L=1:J1
    SK(L)=W*(1-CC(1).^2)/(1+CC(1).^2-2*CC(1)*cos(3.14159*(L-1)/J));
end
if (CC(1) > 0 & CC(1) >= CC(2))
    % 进行白噪声检验
else
    for L=1:J1
        SK(L)=W;
    end
end

% 计算周期、频率、功率谱、95%置信区间
num=1:J+1;
num=num(:)-1;
Period=T(:);
Frequency=1./T(:);
Density=PS(:);
CL95=SK(:);

该函数的主要步骤包括：

1. 计算时间序列信号的自相关系数；
2. 估计粗略的功率谱；
3. 平滑功率谱；
4. 计算功率谱的95%置信区间；
5. 进行红噪声检验，判断信号是否为红噪声；
6. 计算周期、频率、功率谱、95%置信区间。

该函数返回的结果包括周期、频率、功率谱和95%置信区间，可以用于分析时间序列信号的频谱特性。