前回つくった信頼区間つきでスペクトル密度をもとめる関数を、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 件のコメント:
コメントを投稿