ENSO の周期 psd関数のオプション4 平滑化

前回の続き。
これはpsd関数のオプションというわけではないが、Welch法ではなく、得られたスペクトル密度を畳込みによって平滑化してみる。
平滑化にはscipyのcookbookの例にあるsmooth関数を用いる。
平滑化にはスペクトルの周波数の解像度を下げるかわりに統計的自由度を上げる働きがある（見延、2001; Hartmann 2008)。
平滑化には1-2-1フィルタを用いる。numpyでは5点bartlett関数を使うと良い（両端が0）。

import numpy as np
print np.bartlett(5)

[ 0.   0.5  1.   0.5  0. ]

以下がプログラム例。
黒線がもとのスペクトル密度。
赤線が平滑化したスペクトル密度。
（れいによって2年目と7年目の周期に当たるところに縦の点線）

from read_NINO import read_NINO

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


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

#original
Pxx,f=mlab.psd(x, NFFT=len(x),Fs=12)  

# Smooth by 1-2-1 filter 
from cookb_signalsmooth import smooth # from Scipy cookbook
Pxx2=smooth(Pxx[:,0],5,'bartlett') 


#plot
plt.figure()
plt.loglog(f,Pxx,'k-o')
plt.loglog(f,Pxx2,'r-+')
plt.loglog([1./2.,1./2.],[1.E-6,1.E1],'k--') #vertical line for 2 years
plt.loglog([1./7.,1./7.],[1.E-6,1.E1],'k--') #vertical line for 7 years
plt.savefig("NINO_psd_smooth.png")
plt.show()

参考文献
見延　2001:大 気海洋物理学特論4 大気海洋統計データ解析　第3章
Hartmann 2008:ATM552 講義ノート第6a章