スペクトルの推定誤差

見延さんのスペクトル誤差に関するmatlab scriptをpythonに移植してみた。

見延　2001:大 気海洋物理学特論4 大気海洋統計データ解析　第3章
http://ocw.hokudai.ac.jp/Course/GraduateSchool/Science/MeteorologyAndOceanology/2001/page/materials/MeterologyAndOceanology-2001-Note-03.pdf
Page15-16

numpy/scipy/matplotlibを使えば、matlabのscriptを移植するのはそれほど難しくない。
χ二乗分布の累積分布の逆関数はscipyのstats.chi2.ppfをそのまま用いた。
（注：今の場合関係ないが、この関数は自由度が1以下の場合はうまくいかないバグがあるようだ。すくなくともversion 0.7.0で確認した。）

#coding: utf-8
import numpy as np
from numpy.random import randn
import scipy as sp
import scipy.signal as signal

import matplotlib.pyplot as plt 
#---------------------------------------
def chi2inv_nint(p,df):
    import scipy.stats as stats
    c=stats.chi2.ppf(p,df)
    return c

#------------------------------------------
def psd_ex1(tsrs,nfft,ctaper,avfnc):
    """
function [ds_psd,frq]= psd_ex1(tsrs,nfft,ctaper,avfnc)
Power Spectrum Density を直接法によって求めるpython Script
平滑化は周波数領域で行なう．
複数時系列（sz2$gt;1）についても計算する．
入力
tsrs: 入力時系列 (sz1, sz2)各列についてスペクトルを求める
nfft: fft を計算する時間サンプル数
ctaper: taper 関数('hamming','hanning'など)
avfnc: 平滑化ウインドウ
出力
ds_psd: スペクトル密度．１×nfft/2+1．
frq: 周波数．１×nfft/2+1．
todo
PSD の定義に振幅がきちんと合うようにする．
"""
    if tsrs.ndim==1:
        sz1=tsrs.size #tsrs をsz2*1 のデータに変換．
        sz2=1
        tsrs=tsrs.reshape((sz1,1))
    else:
        sz1,sz2=tsrs.shape
    print sz1,nfft
    assert sz1<=nfft, ' sz1$gt;nfft' #nfft はtlng よりも大きくなくてはならない
#満たさないとエラーメッセージを出して終了
    avfnc=avfnc/float(np.sum(avfnc)) #後で説明する平滑化ウインドウの規格化    
    w = getattr(signal, ctaper)(sz1) #実際のtaper 計算
    for n in range(sz2):
        tsrs[:,n]=tsrs[:,n]*w #時系列にテーパーをかける
    tsrs=np.append(tsrs,np.zeros((nfft-sz1+1,sz2)),axis=0) #ゼロパッティング
    famp=sp.fft(tsrs,axis=0) #フーリエ変換
    famp=famp[0:(nfft/2+1),:]/sz1 #実数時系列のスペクトルなら，この範囲のみが必要．
    d_psd=abs(famp)**2; #直接スペクトル値 (d: direct)
    ds_psd=np.zeros_like(d_psd)
    for n in range(sz2):
        ds_psd[:,n]=np.convolve(d_psd[:,n],avfnc,mode='same'); #畳込みを使って平滑化 (ds: direct-smoothed) #畳み込みで増えた部分を除く
    frq=np.arange(nfft/2+1)/float(nfft); #周波数を求めておく
    return ds_psd,frq

########################################################

if __name__ == '__main__':
    plt.figure()
    for npnl in range(4):
        tlng=100; nsmpl=1000
        if npnl==0 :
            nfft=100
            taper=signal.boxcar(tlng)
            ctaper='boxcar'
            avfnc=np.asarray([1, 1, 1, 1, 1]);
            cttl='no padding, boxcar, 5-running'
        elif npnl==1:
            nfft=100
            taper=signal.hamming(tlng)
            ctaper='hamming'
            avfnc=np.asarray([1, 1, 1, 1, 1])
            cttl='no padding, hamming, 5-running'
        elif npnl==2:
            nfft=200
            taper=signal.boxcar(tlng)
            ctaper='boxcar'
            avfnc=np.asarray([1, 1, 1, 1, 1]);
            cttl='double padding, boxcar, 5-running'
        elif npnl==3:
            nfft=200
            taper=signal.hamming(tlng)
            ctaper='hamming'
            avfnc=np.convolve([1, 2, 1],[1, 2, 1])
            cttl='double padding, hamming, double 1-2-1'
    
        tsrs=randn(tlng,nsmpl);
        [ds_psd,frq]=psd_ex1(tsrs,nfft,ctaper,avfnc)
# Monte-Carlo 法で得る90%信頼限界
        srt_psd=np.sort(ds_psd,axis=1)
        c=np.zeros([frq.size,2])                
        c[:,0]=np.log10(srt_psd[:,nsmpl*0.05])
        c[:,1]=np.log10(srt_psd[:,nsmpl*0.95])
        avfnc=avfnc/float(sum(avfnc)) # frequency domain の平滑化関数の規格化
        c0=srt_psd[:,nsmpl*0.5]

        edof=(1+(1/sum(avfnc**2)-1)*tlng/nfft*np.mean(taper))# 等価自由度の計算
        ce=np.zeros([frq.size,2]) 
        ce[:,0]=np.log10(c0*chi2inv_nint(0.050,2*edof)/chi2inv_nint(0.50,2*edof))
        ce[:,1]=np.log10(c0*chi2inv_nint(0.950,2*edof)/chi2inv_nint(0.50,2*edof))
#plot
        plt.subplot(2,2,npnl+1)
        plt.plot(frq,c,'b',frq,ce,'r')
        plt.title(cttl)
        plt.xlabel('frq')
        plt.ylabel('psd')
        if (npnl==0): plt.legend(('Monte-carlo','','Theory'),'lower center',labelspacing=0.05)
    plt.subplots_adjust(wspace=0.6,hspace=0.4)
    plt.savefig("psd_ex1.png")
    plt.show()  

結果の図