2010年4月27日火曜日

NINO3のスペクトル密度を信頼区間つきで



前回つくった信頼区間つきでスペクトル密度をもとめる関数を、NINO3監視指数に適用する。NINO3の読み込みは以前に作ったread_NINOモジュールを用いる。

下の例ではtaperにはHammingフィルタ、周波数空間でのスムージングに1-2-1フィルタ、信頼度は90%を用いた。
黒線がスペクトル密度、赤点線が信頼区間を示す。縦の黒点線は2年と7年に対応する補助線である。


from read_NINO import read_NINO

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import numpy as np

from  psd_withError import psd_withError #import psd_withError module  

data=read_NINO()
x=data["NINO3 ANOM"]

# case1
Pxx,f,upper,lower=psd_withError(x, NFFT=len(x),Fs=12,window=np.hamming(len(x)),smooth=[1,2,1],Confidencelevel=0.90)  
plt.figure()
plt.loglog(f,Pxx,'k-o')
plt.loglog(f,upper,'r--')
plt.loglog(f,lower,'r--')
plt.loglog([1./2.,1./2.],[1.E-5,1.E2],'k--') #vertical line for 2 years
plt.loglog([1./7.,1./7.],[1.E-5,1.E2],'k--') #vertical line for 7 years
plt.savefig("NINO_psd_withError.png")
plt.show() 
 
 

0 件のコメント: