关于平均周期计算的讨论

octopussheng

关于平均周期，本来以为已经很理解了，对于其计算应该也没有任何问题了，但是近期在学习过程中，又发现了很多矛盾的地方，现将我的分析贴出来，也请大家一起帮忙分析分析！

分析对象：Lorenz系统
代码如下：
function dy = Lorenz(t,y,r)
t=16;
r=45.92;
b=4;
dy=zeros(3,1);
dy(1)=-t*(y(1)-y(2));
dy(2)=-y(1)*y(3)+r*y(1)-y(2);
dy(3)=y(1)*y(2)-b*y(3);

求解系统的代码如下：
clear
t0=0;
tf=100;
[t,x]=ode45(@Lorenz,[t0:0.001:tf],[0,1,0]);
结果如图1所示。

在求解功率谱及平均周期时，去掉前面的10002个瞬态项，取后面的稳定结果作为分析的序列，以X序列为对象！
在计算时，用了三种方法，分别是自相关函数法、太阳黑子的例子以及勇俊小论文中的例子。

代码如下：
%power spectrum calculation
%use X as the object
%%%%%%%%%%%%%%%%%%%%%%%%%%%%第一种方法：自相关函数法！
X=x(10002:end,1);
Y=fft(X,1024);
Pyy=Y.*conj(Y)/1024;
f=1000*(0:512)/1024;
plot(f,Pyy(1:513))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%第二种方法：参考自太阳黑子的例子
figure;
YY=fft(X);
N=length(YY);
YY(1)=[];
power=abs(YY(1:N/2)).^2;
nyquist=1/2;
freq=(1:N/2)/(N/2)*nyquist;
plot(freq,power);
period2=1./freq;
figure;
plot(period2,power)
[mp,index]=max(power);
tp=period2(index)
%%%%%最终得到的平均周期：tp =  9000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%第三种方法：参考勇俊：MATLAB在研究非线性混沌中的应用
figure;
YYY=fft(X);
NN=length(YYY);
YYY(1)=[];
power_YYY=log(real(YYY).^2+imag(YYY).^2);
nyquist=1/2;
freq_YYY=(1:NN/2)/(NN/2)*nyquist;
plot(freq_YYY(1:NN/2),power_YYY(1:N/2))
period3=1./freq_YYY;
figure;
plot(period3,power)
[mp,index]=max(power);
tp=period3(index)
%%%%%最终得到的平均周期：tp =  9000
%%%%%%%%%%求平均周期的另一种方法！
%%%sunspot功率谱结果求解平均周期
tp=sum(power)/sum(freq'.*power)
%%结果tp =  970.8212
%%%勇俊例子功率谱结果求解平均周期
tp=sum(power_YYY)/sum(freq_YYY(1:NN/2)'.*power_YYY(1:N/2))
%%结果tp =  9.2646

功率谱图分别如图2～图4所示。
计算平均周期的方法是：太阳黑子中所用方法，即最高功率谱处频率的倒数；刘海龙：基于非线性参数的意识任务分类一文中所用的方法，具体可见图5。
计算结果分别为9000；9000；970.8；9.3，其中自相关函数方法我得到的平均周期居然是inf

这个结果我觉得是不能接受的，也请大家一起参与讨论！解决这个问题！

octopussheng

哎，晕了，看来这个问题没有受到重视啊！！:@L

空山长风

功率谱是自相关函数的傅立叶变换，看你的程序是你先做的傅立叶变换，然后求自相关函数，弄倒了吧

octopussheng

哦，确实是，又查了下书，看来是matlab的帮助里面错了！呵呵，谢谢提醒！