Cross-spectral and power-spectral analysis of real time series, Blackman-Tukey lagged covariance method

実数時系列のクロススペクトル・パワースペクトル解析、 Blackman・Tukey ラグ共分散法

この文書は、増田が執筆した 東京大学大型計算機センター プログラムライブラリ説明書 (東大大型計算機センター, 1987)の原稿をもとに、 操作については、 2000年3月現在増田が使用していたLinuxマシン上での説明に書きかえたものです。 機種やマシンに依存する部分には節見出しに(*)印をつけておきます。 操作以外の部分は1987年の当時の認識に基づいた記述のままです。 ここで採用したクロススペクトルの定義は、 もしかするとスペクトル解析の標準とされているものと (符号のとりかたや規格化の分母などの点で)同じでないかもしれません。 HTMLファイル中の数式およびプログラム例は、<pre>を使っており、 <pre>の中の文字のフォント幅が 日本語文字はASCII文字の2倍でそれぞれ一定であることを前提として書いてあります。 HTML本文中のリンクはこの文書内の相互参照です。

§1. あらまし

1.0 題目

1.1 定義

ここでは定常時系列x(t)の自己共分散関数(ラグτの関数)を

   Cxx(τ) = E(x(t)x(t-τ))                               (1.1)

                ∞
   Sxx(f) = 4 ∫  Cxx(τ)cos 2πfτ dτ                   (1.2)
                0

   Cxy(τ) = E(x(t)y(t-τ))                               (1.3)

   Sxy(f) = Kxy(f) - iQxy(f)；                           (1.4a)

                ∞
   Kxy(f) = 2 ∫ (Cxy(τ)+Cxy(-τ)) cos 2πfτ dτ；     (1.4b)
                0

                ∞
   Qxy(f) = 2 ∫ (Cxy(τ)-Cxy(-τ)) sin 2πfτ dτ       (1.4c)
                0

クロススペクトルの実数部Kxyは、二つの時系列の位相がそろって変化している 振幅を示し、co-spectrumと言われる。 虚数部Qxyは90度位相がずれて変化している振幅を示し、quadrature spectrumと言われる。 位相を問わずに、ある(有限幅の)振動数成分がどの程度同期しているかを論じるには、コヒーレンスという量が使われる。 これはクロススペクトルの絶対値をそれぞれのパワースペクトルの積の平方根で割ったもので、全く無関係のとき0、両者が線形関係になっているとき1、一般にその間の値をとる(このコヒーレンスの2乗をコヒーレンスという人もいる)。 また、二つの時系列の各振動数成分の位相差はクロススペクトルの偏角として求められる(たとえば日野、1977、4章を参照)。

1.2 計算の方式

ここでは、(1.2), (1.4)式に従って、パワースペクトル密度は自己共分散関数の、クロススペクトル密度は相互共分散関数のフーリエ変換として求める。 実際には時系列の長さが有限であるために振動数領域は離散化され、また時系列が離散的にサンプリングされているために振動数領域は有限となる。

特に、共分散を計算するラグ(時間差)LΔtをある有限のところ(L=LMAX)で打ち切ることは、共分散関数に箱型のラグウィンドウ(2.1.3節のIWNTYP=0)をかけていることに相当する。 時間領域でラグウィンドウをかけることと、振動数領域でそのフーリエ変換であるスペクトルウィンドウをたたみこむことは同等であるが、箱型ウィンドウをフーリエ変換すると、非常に波うった形のスペクトルウィンドウとなり、正であるべきパワースペクトルが負になることも多い。 負のパワーを避けるとともにスペクトルの推定誤差を小さくするために、L=LMAXに向かってなめらかに0に落ちるようなラグウィンドウをかけておくことが行なわれる。 ウィンドウを使うことの必要性を誤差解析から体系的に論じたBlackman・Tukey(1958)の本にちなんで、この方法は普通Blackman・Tukey法と呼ばれる。 スペクトルの推定値のばらつきと偏りとは競合する関係にあるので、どんな誤差を最小にしたいかによっていろいろなウィンドウが考えられる。 Hanning (Blackman・Tukeyがオーストリアの気候学者Julius von Hannにちなんでつけた名前)はもともと、(1,2,1)の重みをつけた平滑化をさしていた。 それをスペクトル領域で行なうのと同等なラグウィンドウが2.1.3節のIWNTYP=2である。 スペクトル密度の推定誤差については、正規分布を前提としてではあるが、よく調べられている(たとえば日野、1977、11章参照)。

欠損値(観測もれやデータ伝達の誤りなど)のある時系列のスペクトルの推定については必ずしも統計理論の裏付けはない。 ここでは、共分散関数の推定値を求める際に欠損のある対をはぶくとともに、分母にも算入しないことにする。 つまり、ラグLΔtでの欠損値を含む対の数をmとすると、

                  1     N
   Cxy (LΔt)  = -----  Σ            x(jΔt) y((j-L)Δt)   (1.5)
                 N-L-m  j=L+1(有効な対)

また、共分散関数からスペクトル密度へのフーリエ変換は、

         k         M                          πkL
   Sxx(------) = 4 Σ  Cxx(LΔt) W(LΔt) cos ------  δ(L) Δt
        2MΔt     L=0                          M                (1.6)

         k         M     (e)                  πkL
   Kxy(------) = 4 Σ  Cxy(LΔt) W(LΔt) cos ------  δ(L) Δt
        2MΔt     L=0                          M                (1.7a)

         k         M     (o)                  πkL
   Qxy(------) = 4 Σ  Cxy(LΔt) W(LΔt) sin ------  δ(L) Δt
        2MΔt     L=0                          M                (1.7b)

   δ(L) = 0.5 : L=0とL=M,  
           1.0 : 0 < L < M                                      (1.8)

     (e)      1
   Cxy(τ) = ---(Cxy(τ) + Cyx(τ)): 相互共分散の偶関数部分     (1.9a)
              2

     (o)      1
   Cxy(τ) = ---(Cxy(τ) - Cyx(τ)): 相互共分散の奇関数部分     (1.9b)
              2

Blackman・Tukey法によるスペクトルの推定誤差はほぼ自由度2N/Mのχ2分布に従う。Mを大きくするとスペクトルの分解能1/MΔtは上がるが推定の誤差が大きくなる。 実用的には、MをおよそNの1/10以下にするのがよいとされる(たとえば日野、1977、12章)。

1.3 ほかの計算方式との比較

このほかのスペクトル計算方式としてよく使われるものは、直接フーリエ変換法(普通高速フーリエ変換を使うのでFFT法と呼ばれる)と、最大エントロピー法(MEM)がある。

FFT法は時系列を直接離散フーリエ変換し、振動数領域でかけ合わせることでパワースペクトルやクロススペクトルを求めるものである。 フーリエ変換の一意性から、これは共分散関数から求めるのと同等であることがわかる。 ただし時系列の端の処理だけは互換性がない。 推定誤差を小さくするためのスペクトルウィンドウとしては、Blackman・Tukey法で考えられているラグウィンドウをフーリエ変換したものを使えばよい。 ほしい振動数の分解能が高くなるにしたがってBlackman・Tukey法ではラグを長くとらなければならないので共分散の計算量がふえ、相対的にFFT法が有利になる。ただしFFT法では欠損値の処理が困難である。

MEM法あるいは自己回帰モデルあてはめ法は、時系列が線スペクトルでよく近似されるような周期性をもっていることを前提として、それを情報量を失わないように取り出そうとするものである。 それに対してBlackman・Tukey法やFFT法は、仮定をおかずにある振動数帯にどれだけのパワーをもっているかを見ようとするものである。 したがって、MEMに比べてBlackman・ Tukey法ではするどいピークが出にくいのは当然である。

1.4 ルーチンの構成

G3RIL  G3RXLまたはG3RPLに必要な三角関数とラグウィンドウを用意する。
G3RXL  二つの時系列のクロススペクトルを計算する。
G3RPL  一つの時系列のパワースペクトルを計算する。
G3RAV  平均値を計算する(G3RXL,G3RPLの下うけ)。
G3RCV  共分散を計算する(G3RXL,G3RPLの下うけ)。
G3RFU  共分散からスペクトル密度へのフーリエ変換(G3RXL,G3RPLの下うけ)。
Y8AA01 二つの時系列のクロススペクトルを計算し、印刷する(G3RXLの応用)。

  Y8AA01----G3RXL----G3RAV
                 ----G3RCV
                 ----G3RFU
        ----G3RIL
            G3RPL----G3RAV
                 ----G3RCV
                 ----G3RFU

§2.使い方

2.0 コマンド(*)

2.0.1 ロードモジュールライブラリの利用

% f77 USER.f -lg3rxl
% ./a.out

USER.fの内容

      PROGRAM p1
..........
      CALL G3RIL(......)
      CALL G3RXL(......)
..........
      END

注: ルーチンY8AA01はロードモジュールライブラリに含めていない。

2.0.2 ソースライブラリの利用(*)

2.0.2.1 プログラムの実行(*)

% cp /usr/local/src/masuda/spectr/g3rxl.f .
% f77 USER.f
% ./a.out

USER.fの内容

      PROGRAM p2
..........
      CALL G3RIL(......)
      CALL G3RXL(......)
..........
      END
      INCLUDE '/usr/local/src/masuda/spectr/g3rxl.f'

注: ルーチンY8AA01を利用する場合は、 g3rxl.fに加えて、同じディレクトリの別ファイルy8aa01.fも INCLUDEする必要がある。

2.0.2.2 ライブラリプログラムの印刷(*)

cd /usr/local/src/masuda/spectr
a2ps g3rxl.f | lpr -Pgslips

2.0.2.3 ライブラリプログラムの複写(*)

cp /usr/local/src/masuda/spectr/g3rxl.f .

2.1 初期設定ルーチン G3RIL

2.1.0 役割

2.1.1 呼びだし方

CALL G3RIL(LMAX,MXCS,COSSIN,WINDOW,IWNTYP,LDATA,QLAMDA)

2.1.2 引き数の説明

2.1.3 ラグウィンドウの選択 IWNTYPの意味

2.1.4 ラグウィンドウの補正 LDATAの意味

                  1    N
   Cxy (LΔt) = -----  Σ x(jΔt) y((j-L)Δt)                    (2.1)
                 N-L  j=L+1

                  1    N
   Cxy (LΔt) = ----- Σ x(jΔt) y((j-L)Δt)                    (2.2)
                  N   j=L+1

2.2 クロススペクトルの計算 G3RXL

2.2.0 役割

2.2.1 呼び出し方

CALL G3RXL(LENG,  LMAX,  DELT,  A1,    A2,
     COV1,  COV2,  COV12, COV21, COVC,  COVS,  MISS1, MISS2, MIS12, MIS21, 
     FOUR1, FOUR2, FOURC, FOURQ, FOURX, FOURY, COHER, PHASE,
     MXCS,  COSSIN,WINDOW,IFMIS1,QMISS1,IFMIS2,QMISS2,IRTN)

2.2.2 引き数の説明

2.2.3 エラーメッセージ

  G3RXL: SPECTRA ARE NOT CALCULATED.
  ALL DATA ARE MISSING FOR COVARIANCE WITH LAG = (ラグの数)
  BETWEEN SERIES NO.(1または2)AND(1または2)

2.3 パワースペクトルの計算 G3RPL

2.3.0 役割

2.3.1 呼び出し方

CALL G3RPL(LENG,  LMAX,  DELT,  A1,    COV1,  MISS1,
           FOUR1, MXCS,  COSSIN,WINDOW,IFMIS1,QMISS1,IRTN)

2.3.2 引き数の説明

2.3.3 エラーメッセージ

  G3RPL: SPECTRA ARE NOT CALCULATED.
  ALL DATA ARE MISSING FOR COVARIANCE WITH LAG = (ラグの数)

2.4.1 呼び出し方

2.4.2 引き数の説明

2.5 共分散関数の計算 G3RCV

2.5.0 役割

2.5.1 呼び出し方

CALL G3RCV(LENG,LMAX,A1,A2,COVX,MISSX,ABAR1,ABAR2,
           IFMIS1,QMISS1,IFMIS2,QMISS2,IRTN)

2.5.2 引き数の説明

2.6 共分散関数からスペクトルへの変換 G3RFU

2.6.0 役割

2.6.1 呼び出し方

2.6.2 引き数の説明

2.7 クロススペクトルの計算と印刷 Y8AA01

2.7.0 役割

2.7.1 呼び出し方

2.7.2 引き数の説明

2.7.3 印刷されるもの

1. 三角関数とラグウィンドウの値 (KAISU = 0を与えた場合):
   それぞれのKLの値ごとに KL, COSSIN(KL,1), COSSIN(KL,2),
   それぞれのLの値ごとに (実際には上の三角関数の横に) L, WINDOW(L) を表示する。
2. 共分散関数:
   それぞれのLの値ごとに L, DELT*L, COV1(L), MISS1(L), COV2(L), MISS2(L), COV21(L), MIS21(L), COV12(L), MIS12(L), COVC(L), COVS(L), L を表示する。
3. クロススペクトル:
   それぞれのKの値ごとに K, 振動数FREQ=K/(2*LMAX*DELT), 周期PERIOD=1/FREQ, FOUR1(K), FOUR2(K), FOURX(K), FOURC(K), FOURQ(K), FOURY(K), COHER(K), 度単位の位相差PHASE(K)*180/π,周期単位の 位相差PHASE(K)/2π, K を表示する。

2.7.4 制限

   Y8AA01: LMAX (LMAXの値) IS LARGER THAN DIMENSION (MXLの値)

   PARAMETER (MXL=100, MXCS=2*MXL-1)

2.7.5 Common変数(*)

   COMMON   COSSIN(0:MXCS,2), WINDOW(0:MXL)

2.7.6 仮定されている値

IWNTYP = 2           Hanningである。
LDATA  = 0           2.1.4節で述べた補正を行なわない。
QLAMDA = 0.0         (Hanningなので特に意味はない。)

IFMIS1 = IFMIS2 = 1      どちらの時系列にも欠損値がありうる。
QMISS1 = QMISS2 = 1.E6   欠損値は絶対値百万以上の値で表わされている。

2.7.7 隠しパラメータ KSKIP

サブルーチンY8AA01の中にKSKIPという変数がある。登録されているプログラムでは KSKIPはDATA文で1に設定してあるが、これを0に変えてコンパイル、実行すると、2.7.3節で述べたものに加えて、クロススペクトルを振動数の小さいほうから積分したものおよびそれを全パワーで規格化したものを表示する。 なお、旧版では、規格化の際振動数の番号0、1、2の部分を別あつかいにしていた(おそらく低周波部分はスペクトルの信頼度が落ちるため)が、一般性がないので、全パワーに対する割合に改めた。

§3. 使用例

3.1 tg3rpl.f (*)

% cp /usr/local/src/masuda/spectr/tg3rpl.f .
% f77 tg3rpl.f -lg3rxl
% ./a.out

3.2 tg3rxl.f (*)

% cp /usr/local/src/masuda/spectr/tg3rxl.f .
% f77 tg3rxl.f -lg3rxl
% ./a.out

3.3 ty8aa01.f (*)

% cp /usr/local/src/masuda/spectr/ty8aa01.f .
% cp /usr/local/src/masuda/spectr/y8aa01.f .
% f77 ty8aa01.f y8aa01.f -lg3rxl
% ./a.out

§4. 備考

4.1 経過

このプログラムの原型は、1960年代後半に東京大学理学部地球物理学教室の気象研究室で、熱帯大気の波動擾乱のデータ解析に使われていたサブルーチンSPECTRである。 直接にはMunkら(1959)の海洋の大規模波動のデータ解析の論文とPanofsky・Brier (1958)の気象統計学の教科書を参考にしているが、いずれもスペクトル解析に関してはTukey(1949)に従っており、Blackman・Tukey (1958)の本で理論的裏付けが与えられたのと同じ方法である。

このプログラムはMaruyama(1968)をはじめとする多くの研究に使われた。 遠藤(1971)はこれをソースリストつきで紹介している。 また、東大大型計算機センターのプログラムライブラリにもY8/TC/AA01として1972年に登録された(東大大型計算機センター、1973)。 [作成: 林 良一(Yoshikazu HAYASHI)、新田 勍(Tsuyoshi NITTA), 1972- 3-31]

当時はデータはカードで入力し、結果はラインプリンタに出てくるのが普通だったが、その後計算機の使い方が変わり、データはおもに磁気テープの形で提供され、TSS(会話型処理)で計算し、結果はグラフィック装置に作図するのが普通になった。 SPECTRもプリンタ出力を前提としたそのままの形では使われず、結果を引き数の形で持ち帰るように改造して使うことが多くなった。

そこで、1986年に次の改訂を行なった。 [改訂: 増田耕一 (Kooiti MASUDA), 1987- 1- 6]. SPECTRのうち実際に計算を行なっている部分をG3RXLとして分離し、利用者が直接呼べるようにした。 単独の時系列のパワースペクトルを計算する需要も大きいので、そのためのルーチンG3RPLも用意した。 初期設定ルーチンG3RILを整備し、Hanning以外のいろいろなラグウィンドウが選択できるようにした。 全体に、Fortran77の構文を使ってプログラムの書体を読みやすいものにするよう努力した。 サブルーチンY8AA01は従来のSPECTRの形で使えるように残したので、以前からの利用者はルーチン名とコマンド名だけ変えればよい。

4.2 参考文献

Akaike, Hirotugu, 1962:

On the design of lag window for the estimation of spectra. Annals of the Institute for Statistical Mathematics, Vol. 14, 1 - 21.

赤池弘次,1968:

スペクトル解析。 磯部 孝(編), 相関函数およびスペクトル, 東京大学出版会、28 - 46.

赤池弘次,中川東一郎,1972:

ダイナミックシステムの統計的解析と制御。サイエンス社。 189 pp.

Blackman, R. B. and Tukey, J. W., 1958:

The Measurement of Power Spectra. Dover. 190 pp.

遠藤昌宏, 1971:

気象擾乱のスペクトル解析。気象研究ノート(日本気象学会), 110号、134 - 145.

日野幹雄, 1977:

スペクトル解析。 朝倉書店。 300 pp.

藤井光昭, 1974:

時系列解析。 コロナ社。222 pp.

Maruyama, Taketo, 1968:

Time sequence of power spectra of disturbances in the equatorial lower stratosphere in relation to the quasi-biennial oscillation. Journal of the Meteorological Society of Japan, Vol. 46, 327 - 342.

Munk, W. H., Snodgrass, F. E. and Tucker, M. J., 1959:

Spectra of low-frequency ocean waves. Bulletin of Scripps Institution of Oceanography, Vol. 7, 283 - 362.

Panofsky, H. A. and Brier, G. W., 1958:

Some Applications of Statistics to Meteorology. Pennsylvania State University, 224 pp.

東京大学大型計算機センター, 1973:

Y8/TC/AA01, 擾乱のパワースペクトル解析。 東京大学 大型計算機センターニュース、Vol. 8, Supplement 2, 246 - 251.

東京大学大型計算機センター, 1987:

G3/TC/RXL, 実数時系列のクロススペクトル・パワースペクトル解析、Blackman・Tukey ラグ共分散法。 東京大学 大型計算機センターニュース、Vol. 19, Supplement 1, 69 - 85.

Tukey, J. W., 1949:

The sampling theory of power spectrum estimates. Symposium on Applications of Autocorrelation Analysis to Physical Problems, Woods Hole, June 13, 1949. NAVEXOS-P-735, Office of Naval Research.