如何讨论一个时间序列的功率谱

qqqin

在这篇文献中有这样一张图（原文图3），如下所示，其中这个非线性系统由几个子系统（lorenz63与slab upper ocean ，simple pycnocline predictive model）耦合而成。通过利用matlab的ode45可以产生对应每个系统状态量的时间序列数据， 其中，图片所示的其中3个状态量的功率谱(x_2,omega,eta).按我的理解，这个功率谱可通过先对时间序列求自相关，然后对自相关进行fft求得，但是结果看起来不太一样。
不知道这样的处理过程是否有错误，或者哪里有遗漏？同时，也不太理解其中的95% cconfidence-upper-limit是怎么画出来的？还望了解相关知识的指点一二，谢谢了.

1. sigma=9.95;k=28;b=8/3;C1=0.1;C2=1;Od=1;Om=10;Sm=10;
   
2. Ss=1;Spd=10;Sigma=100;C3=0.01;C4=0.01;C5=1;C6=0.001;
   
3. ntime=4000;
   
4. dt=0.001; T=dt*ntime;
   
5. initial=[0,1,0,0,0]';
   
6. tspan=dt:dt:dt*ntime;
   
7. [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);
   
8. %states(:,1),states(:,2),states(:,3),states(:,4),states(:,5)  == x_1,x_2,x_3,omega,eta
   
9. %need some code to reproduce the power spectrum of Fig3
   
10. %and the confidence upper limit
    
11. %consider fft
    
12. function dy = couple(y,t,sigma,k,b,C1,C2,Od,Om,Sm,Ss,Spd,Sigma,C3,C4,C5,C6)
    
13.     dy(1)=-sigma*y(1)+sigma*y(2);
    
14.     dy(2)=-y(1)*y(3)+(1+C1*y(4))*k*y(1)-y(2);
    
15.     dy(3)=y(1)*y(2)-b*y(3);
    
16.     dy(4)=(C2*y(2)+C3*y(5)+C4*y(4)*y(5)-Od*y(4)+Sm+Ss*cos(2*pi*t/Spd))/Om;
    
17.     dy(5)=(C5*y(4)+C6*y(4)*y(5)-Od*y(5))/Sigma;
    
18.     dy=[dy(1);dy(2);dy(3);dy(4);dy(5)];
    
19. end
    

复制代码

something210

我也看过这篇文献，没有重复出来····

qqqin

something210 发表于 2019-1-4 09:29
我也看过这篇文献，没有重复出来····

你是怎么处理的呢，一起讨论下