2010年3月12日金曜日

ENSOの周期



ENSO(エルニーニョ)の周期を見てみる。

データは前回作ったモジュールでNINO3領域ののアノマリを得る。

スペクトル密度を求める関数がmatplotlib.mlab(Matlab風の関数群)にpsdとしてある。
(図まで描いてくれるmatplotlib.pyplotのpsdというのもある。)

そのpsdを用いて作ったスペクトル密度の図とプログラムが以下の通り。
データ数(NFFT)はすべてを使う。
時間の単位は年とする。毎月のデータがあるから単位時間あたりのサンプル数(Fs)は12。
エルニーニョは2年から7年の周期があると言われているから補助線として2年(周波数1/2)と7年(周波数1/7)に当たるところに補助線(点線)を引いている。実際、それらの間にスペクトルピークが見られている。
from read_NINO import read_NINO

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt

data=read_NINO()
x=data["NINO3 ANOM"]
Pxx,f=mlab.psd(x, NFFT=len(x),Fs=12)  
plt.figure()
plt.loglog(f,Pxx)
plt.loglog([1./7.,1./7.],[1.E-9,1.E2],'k--')
plt.loglog([1./2.,1./2.],[1.E-9,1.E2],'k--')
plt.savefig("NINO_psd.png")
plt.show()

0 件のコメント: