О спектре мощности и верхнем пределе достоверности данных временного ряда

Question

Пока у меня есть связанная система с 5 переменными, и я использую метод Рунге-Кутты для интеграции.

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=4000;
%X=zeros(5,ntime);
dt=0.001; T=dt*ntime;
%X(:,1)=[0,1,0,0,0]';
initial=[0,1,0,0,0]';
tspan=dt:dt:dt*ntime;
%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
[time,states] = ode45(@(t,y) couple(y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6),tspan, initial);


 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;
    dy=[dy(1);dy(2);dy(3);dy(4);dy(5)];

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+dt/2,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
    k2 = dt*dy;
    dy = couple(y+k2/2,t+dt/2,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6);
    k3 = dt*dy;
    dy = couple(y+k3,t+dt,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

И я пытаюсь использовать fft рассчитать спектр мощности

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

cxn = xcorr(xn,'unbiased'); 
CXk = fft(cxn);
psd2 = abs(CXk);
psd2 = psd2/max(psd2);
psd2 = 10*log10(psd2);
% index = 0:round(nfft/2-1);
% k = index*Fs/nfft;
% psd2 = psd2/max(psd2);
% psd2 = 10*log10(psd2(index+1));
% figure(1)
% plot(k,psd2);
plot(psd2)
xlim([0 100])
grid on

но результаты кажутся различными с помощью ссылки ( https://www.nonlin-processes-geophys.net/24/681/2017/npg-24-681-2017.pdf). С другой стороны, я не совсем понимаю, как получить верхний предел достоверности, как показано на рисунке 3. Надеюсь, что кто-нибудь может дать мне несколько советов. Заранее спасибо.