ENSOの周期 psd関数のオプション3 Welch法

前回の続き。
ここで用いているpsd関数は、データ長をブロックにわけ、その平均をおこなうことでスペクトル密度をスムーズにするオプションがある(Welch法）。

The power spectral density by Welch’s average periodogram method. The vector x is divided into NFFT length blocks.
noverlap: integer The number of points of overlap between blocks. The default value is 0 (no overlap).

以下がその例。
黒線がオリジナル。
赤線が、データ長を二分割したケース（オーバラップ0）。
青線が、データ超を二分割した上で、50%のオーバラップを許した例。
スペルトル密度の線がなめらかになるのがわかる。
(補助線として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)  

# Date are devided into 2 blocks
#Zero pad to the length of the orignal data sequence.
Pxx2,f2=mlab.psd(x, NFFT=len(x)//2,pad_to=len(x),Fs=12)

#Date are devided into 2 blocks
#Allow 50 % overlap 
#Zero pad to the length of the orignal data sequence.
Pxx3,f3=mlab.psd(x, NFFT=len(x)//2,noverlap=int(0.5*len(x)/2.),pad_to=len(x),Fs=
12)

#plot
plt.figure()
plt.loglog(f,Pxx,'k-o')
plt.loglog(f2,Pxx2,'r-+')
plt.loglog(f3,Pxx3,'b-x')
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_welch.png")
plt.show()

参考
pylab.psd関数のデモ
pylab_examples example code: psd_demo2.py
MATLABのWelch法の関数pwelch
http://www.mathworks.com/access/helpdesk_ja_JP/help/toolbox/signal/pwelch.html