关于时间序列的功率谱
功率谱表示随机运动过程在各频率成分上的统计特性，是研究随机振动的基本工具。对于给定的随机信号，可以采用标准程序软件计算货专用的频谱分析仪测定其功率谱。为了描述混沌振动的随机性，可以应用研究随机振动的频谱分析方法识别混沌振动。通常假设混沌是各态历经的，即时间上的平均量与空间上的平均量相等。
在试验测量和计算机仿真中，得到的往往是相差相同时间间隔τ的时间序列：
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=;
求解系统微分方程程序：
clear
ff=0.222;
options=odeset('RelTol',1e-7);
t0=0;
tf=100;
=ode45(@dafen,,,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;
=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));
结果图比较如下：
虽然功率谱这个东东论坛里面讨论的也很多了，我还是迎大家一起参与讨论，把这个东东搞透彻！
附件:
http://forum.vibunion.com/data/attachment/album/2007/09/51118_200709281550391.thumb.jpg
duffing系统的相图[时间:2007-9-28 15:50]
http://forum.vibunion.com/data/attachment/album/2007/09/51118_200709281550591.thumb.jpg
x(1)的时间历程[时间:2007-9-28 15:50]
http://forum.vibunion.com/data/attachment/album/2007/09/51118_200709281551221.thumb.jpg
方法一得到的功率谱图[时间:2007-9-28 15:51]
http://forum.vibunion.com/data/attachment/album/2007/09/51118_200709281551411.thumb.jpg
方法二采用periodogram函数的结果图[时间:2007-9-28 15:51]
http://forum.vibunion.com/data/attachment/album/2007/09/51118_200709281552051.thumb.jpg
方法二采用fft后求平方和的结果图[时间:2007-9-28 15:52] 你这个时间序列是在用数值解的过程中，时间间隔是步长吗？不是数据采集过程中的采样时间间隔吗？我还没有转过来，计算过程中设计的时间序列，有实际意义吗？求出它的功率普有什莫用呢？ 请问还有别的方法吗,各种方法都针对哪类系统效率更好?
回复 #2 vib 的帖子
采样时间间隔和时间步长是不一样的对于你后面的问题我不是很理解，能不能解释一下哈！
回复 #2 vib 的帖子
数值仿真的结果一般是时间历程，所以FFT分析是想知道它的一些频域特征，所以计算功率谱回复 #3 skywm 的帖子
关于功率谱的计算基本上就是上面的几种方法了，如果还有其他的更好的方法，也欢迎提出来一起分享！ 想问下楼主 作出的功率谱纵坐标单位是什么呢?和伯德图有没有什么相似的地方? 纵坐标的单位和你想做的如“位移”、“速度”、“加速度”等的不同而不同，关于功率谱密度的单位在振动板块有一个帖子的，你可以搜索参考一下！ 楼主在程序中讲要 略去前面的一些瞬态解，并且这个数量可以自己按照求解精度来确定，如何确定呢，能不能说的具体点啊？谢谢 请问楼主,xx1是什么数据? 请问你的时间序列是几维的？如果是四维的，怎么求功率谱? 程序中的nfft的值如何确定？ nyquist=1/2;freq1=(1:N1/2)/(N1/2)*nyquist;
这句程序就决定了频率的最大值是0.5，是不是不对啊？
页:
[1]