ENSO の周期 psd関数のオプション5 scale_by_freq

前回の続き。
これは単なるテクニカルなオプションだが、メモのために記しておく。

*scale_by_freq*: boolean
          Specifies whether the resulting density values should be scaled by the scaling frequency, which gives density in units of Hz^-1. This allows for integration over the returned frequency values. The default is True for MatLab compatibility.

今まで使ってきた下のような例では、もとは月単位のデータを年単位のデータにしたときに(Fs=12)、自動的に12で割っている（scale_by_freq=Trueがデフォルト)。
すなわち、周波数が(1/2,1/3,1/4, ......)（月単位）の時のスペクトル密度（Pxx2,scale_by_freq=False)の場合と、周波数が（6,4,3, ....)(年単位）の時のスペクトル密度Pxx(scale_by_freq=True)の積分が一致するように

Pxx=Pxx2/12

の関係がある。
下のグラフのようにPxx(黒線）とPxx2/12(赤線）とは一致する。

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"]
Pxx,f=mlab.psd(x, NFFT=len(x),Fs=12)  
Pxx2,f2=mlab.psd(x, NFFT=len(x),Fs=12,scale_by_freq=False) 

plt.figure()
plt.loglog(f,Pxx,'k')
plt.loglog(f2,Pxx2/12.,'r--')
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_scale.png")
plt.show()