功率谱表示随机运动过程在各频率成分上的统计特性，是研究随机振动的基本工具。对于给定的随机信号，可以采用标准程序软件计算货专用的频谱分析仪测定其功率谱。
为了描述混沌振动的随机性，可以应用研究随机振动的频谱分析方法识别混沌振动。通常假设混沌是各态历经的，即时间上的平均量与空间上的平均量相等。

在试验测量和计算机仿真中，得到的往往是相差相同时间间隔τ的时间序列：
x(1),x(2),...,x(N)
方法一：
对该序列附加周期性条件x(N+i)=x(i)(i=1,2,...)后可以计算相关函数c(i)，在matlab中采用xcrr命令可直接求解，然后对结果作离散傅立叶变换，得到的结果p(j)即表示x(k)中第j个频率成分，即为时间序列的功率谱。
方法二：
不经求解自相关函数，而直接求解x(i)的离散傅立叶系数a(j)和b(j)（计算公式这里就不给出了，很多书上都有的）,然后计算
p(j)＝a(j)^2+b(j)^2)
通常为许多组{x(i)}计算一批{p(j)}，平均后即逼近序列的功率谱。

下面以duffing系统为例进行求解说明！
系统：x''=f*cos(1.2*t)-x^3+x-0.3*x'

定义系统微分方程程序：
function df=dafen(t,x,flag,force)
df=[x(2);force*cos(1.2*t)-x(1)^3+x(1)-0.3*x(2)];

求解系统微分方程程序：
clear
ff=0.222;
options=odeset('RelTol',1e-7);
t0=0;
tf=100;
[t,x]=ode45(@dafen,[t0:0.001:tf],[0,0],options,[],ff);
plot(x(10002:end,1),x(10002:end,2))

这里对时间序列的获取进行一些说明：
求解微分方程的时候，步长h=0.001,在略去前10002个瞬态解（这个数量可以自己按照求解精度来确定）后，取采样时间τ=0.1s,共选取900个数据，即序列长度N=900，显然采样频率f=100。

方法一的代码：
%第一种方法：
%先对时间序列求自相关函数
%然后对自相关函数进行傅立叶变换，即可得到时间序列的功率谱
Fs=100;%采样频率
N=length(xx1);
nfft=1024;
cxx1=xcorr(xx1);%计算序列的自相关函数
CXX1k=fft(cxx1,nfft);%求功率谱
Pxx1=abs(CXX1k);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx1=10*log10(Pxx1(index+1));
plot(k,plot_Pxx1)
同样，x(2)的功率谱也可以绘制了！

方法二的代码：
采用方法二的时候，我比较了一下用periodogram函数和直接求平方和的方法之间的区别，代码分别如下：
采样窗函数
Fs=100;%采样频率
window1=boxcar(length(xx1));%矩形窗
nfft=1024;
[Pxx1,f1]=periodogram(xx1,window1,nfft,Fs);
plot(f1,10*log10(Pxx1))

fft之后求平方和：
y1=fft(xx1);
N1=length(y1);
y1(1)=[];
power1=log(real(y1).^2+imag(y1).^2);
nyquist=1/2;
freq1=(1:N1/2)/(N1/2)*nyquist;
plot(freq1(1:N1/2),power1(1:N1/2));

结果图比较如下：
虽然功率谱这个东东论坛里面讨论的也很多了，我还是迎大家一起参与讨论，把这个东东搞透彻！