ARモデルのPSD計算法
自己回帰モデル（AR法）を用いたパワースペクトル密度算出　その３

ここでは、（１）式の自己回帰モデル（ARモデル）の次数Mやパラメータが決定されたとして、モデルパラメータを使って（２）式のようにパワースペクトル密度が計算されることを解説します。

（１）

（２）

ARモデルのパラメータをユールウォーカ法（最尤法）で推定する方法は「ARモデルのパラメータ推定法」を、ARモデルの骨格の選択（ARモデルの次数選択）は「AICによるARモデルの次数決定法」を参考にしてください。

ARモデルのインパルス応答表記

はじめに、時間遅延演算子（ラグ演算子ともいいます）を導入します。この演算子は、時系列データを「インパルス応答」とよばれる形式で記述するために必要となってくるものです。ここまでに、時系列データをARモデルで記述する手順を解説しましたが、そのARモデルをさらにインパルス応答の形式で記述してみよう、ということです。

時間遅延演算子Zは（３）で示されるように、時間を１単位づつ遅らせる演算子です。

（３）

時点ｓのｘの値に対して、Zをｊ回演算すると、時点ｓ－ｊのｘの値が得られる、というのが（３）式の示していることです。この演算子Zを使って（１）式のARモデルを書き直すと式（４）になります。

（４）

Zは演算子でしたが、普通の変数と同じように扱えるものとして、（４）式をについて解きますと（５）式になります。変数ｘに作用する演算子として導入された時間遅延演算子Zを、普通の変数と同様の扱いをしてよいものか疑問に思う方もあると思いますし、それは正しい疑問です。が、ここでは、演算子を変数のように扱うと便利だし上手くいくからやってみた、程度に軽く考えて使い方に慣れてください。

（５）

（５）式から、式変形（多項式の割り算）すると、（６）式のように係数ｈの級数に変換できます。

（６）

Zは時間遅延演算子として同様、にも作用し、のような作用になります。(５)式から（６）式への式変形は多項式の割り算によって、係数ｈが以下のように計算できる事が確認できます。ただし、です。

簡単な例で確かめてみます。Zの時数Mは２とします。

実際に計算してみると、係数ｈは確かに以下のように表記できることがわかる。

補足ですが、（６）式の結果を時間遅延演算子Zを使わないで表記すると、（７）式です。

（７）

以上のように、式（１）で表現されるARモデルの時系列データは、（６）式、（７）のように時刻ｓまでに加わったノイズ（は分散のホワイトノイズ）の重ね合わせで表現できることがわかった。この（６）式、（７）式の表現形式を「インパルス応答」と呼ばれ、デジタル信号処理技術の分野でよくつかわる形式です。

 

ARモデルのパワースペクトル

インパルス応答をつかったPSD表記（ウィーナーヒンチンからの導出）

「パワースペクトル密度と自律神経バランス指標LF/HF」の（９）式に示した、自己相関関数のフーリエ変換がパワースペクトルになるという、ウィナーヒンチンの定理がスタートとなります。上記では、連続時間の場合の結果でしたが、離散時間の場合のウィーナーヒンチン定理は（８）式になります。

（８）

自己相関関数C(ｋ)は、その定義から、時間がｋだけ離れたｘの信号同士の掛け算の期待値ですから、となります。ここで、時系列データｘはホワイトノイズ系列を入力としたインパルス応答として記述できるという、先の（７）式の結果を利用して（８）式を式変形します。

（９）

（９）式は和の記号が並ぶ複雑な形式になってしまいましたが、ホワイトノイズの時系列であるの次の特徴を利用することを考えます。同士の積の期待値は、同じ時刻ｓどうしの積であれば、違う時刻どうしの積であれば０です。

（１０）

これに注意すると、（９）式の　　の部分は j = k+l の部分のみが値を持ち、その他組み合わせの時は０なので簡単化できます。

このように非常に単純な形になります。さらに周期関数の部分を２つに分解して、整理しなおすと、

（１１）

パワースペクトルは、インパルス応答を使って１１式のように書けることがわかりました。

（１１）式の　　の部分はよく見ると非常に直観的にわかりやすい形になっています。つまり、インパルス応答（時系列データはこれにホワイトノイズが入力であるときの出力でした。）に周期fの周期関数を入力したときの出力がどのくらい増幅、あるいは減衰するかを示しています。この部分を周波数応答関数と呼びます。

インパルス応答とARモデル係数の関係からARモデルのパワースペクトル密度を導出

ここでもう一度、（６）式の結果を使います。時間遅延演算子Zを普通の変数と同様に扱って、多項式の割り算をした結果導かれた関係式でした。必要な部分だけを抽出して書くと（１２）式になります。

（１２）

（１２）式の関係は、Zにどのような値を入れても成り立つはずで、　を代入すると、

（１３）

（１３）式のようにARモデルの記述へと変換することができます。実用上も、ARモデル係数が求まっているのであればインパルス応答形式を経由せずに、直接（１３）式を計算したほうが簡単に済みます。以上で、ARモデルの係数がわかれば直ちにパワースペクトルが（13)式で求まることが確認できました。

matlab (Signal Processing Toolbox が必要)　には、時系列データからユールウォーカ法でARモデルを求めてパワースペクトルを計算するpyulear(vector,order) が実装されています。具体例はこちらに掲載したサンプルコードを参考にしてください。