# About the power spectrum and confidence upper limit

For now, I have a coupled system with 5 variables and use the Runge-Kutta method to integrate.

sigma=9.95;k=28;b=8/3;C1=0.1;C2=1;Od=1;Om=10;Sm=10;
Ss=1;Spd=10;Sigma=100;C3=0.01;C4=0.01;C5=1;C6=0.001;

ntime=1000;
X=zeros(5,ntime);
dt=0.001; T=dt*ntime;
X(:,1)=[0,1,0,0,0]';
for t = 2:ntime
y = X(:,t-1);
X(:,t)=rk4(dt,y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
end

function dy = couple(y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6)
dy(1)=-sigma*y(1)+sigma*y(2);
dy(2)=-y(1)*y(3)+(1+C1*y(4))*k*y(1)-y(2);
dy(3)=y(1)*y(2)-b*y(3);
dy(4)=(C2*y(2)+C3*y(5)+C4*y(4)*y(5)-Od*y(4)+Sm+Ss*cos(2*pi*t/Spd))/Om;
dy(5)=(C5*y(4)+C6*y(4)*y(5)-Od*y(5))/Sigma;

end

function y_next = rk4(dt,y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6)
dy = couple(y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
k1 = dt*dy';
dy = couple(y+k1/2,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
k2 = dt*dy';
dy = couple(y+k2/2,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
k3 = dt*dy';
dy = couple(y+k3,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
k4 = dt*dy';
y_next = y+k1/6+k2/3+k3/3+k4/6;
end


And I try to use fft to calculate the power spectrum

xn = X(4,:); % X(1,:)X(2,:)X(3,:) for x_1,x_2,x_3, X(4,:) for omega, X(5,:) for eta.
Fs= 1/dt;
nfft = 2 ^ nextpow2(ntime);
cxn = xcorr(xn,'unbiased');

CXk = fft(cxn);
psd2 = abs(CXk);

index = 0:round(nfft/2-1);
k = index*Fs/nfft;
psd2 = psd2/max(psd2);
psd2 = 10*log10(psd2(index+1));
plot(k,psd2);
xlim([0 100])
grid on


but the results seem different with the reference(https://www.nonlin-processes-geophys.net/24/681/2017/npg-24-681-2017.pdf). On the other hand I don't quite understand how to acquire the confidence upper limit as the reference Figure3 shown. Hope that anyone can give me some advice. Thanks in advance.