mem.py 技術ドキュメント

プログラムの動作

mem.py は、時系列データに対して最大エントロピー法（Maximum Entropy Method, MEM）を適用し、その周波数スペクトルを解析するためのPythonプログラムです。特に、Burg法を用いて自己回帰（AR）モデルの係数を推定し、それに基づいてパワースペクトルを計算します。

主な機能は以下の通りです。

  1. 時系列データ生成: 複数のサイン波とランダムノイズを合成して、解析対象となるサンプル時系列データを内部で生成します。

  2. 平均値除去: 解析前に時系列データから平均値を除去し、直流成分の影響を排除します。

  3. MEMによるスペクトル解析: Burg法とLevinson-Durbin再帰アルゴリズムを用いて、入力時系列データのパワースペクトル、予測誤差パワー、および自己相関関数を計算します。

  4. 結果の出力: 計算された周波数、パワースペクトル、予測誤差パワー、自己相関関数を標準出力に表形式で表示します。

  5. グラフ描画: 入力時系列データと計算されたパワースペクトルを matplotlib を用いてグラフ化し、画面に表示します。

このプログラムは、時系列データの周期性や隠れた周波数成分を高分解能で検出したい場合に有用です。

原理

mem.py は、最大エントロピー法（MEM）に基づいて時系列データのパワースペクトルを推定します。MEMは、限られた自己相関情報から、それと整合しつつエントロピーが最大となるようなパワースペクトルを推定する手法です。これにより、フーリエ変換などの古典的な方法と比較して、短いデータ長でも高い周波数分解能が得られるという利点があります。

本プログラムでは、MEMの係数推定にBurg法を採用しています。Burg法は、前向きおよび後ろ向き予測誤差の二乗和を最小化することで、自己回帰（AR）モデルの係数（パーシャル自己相関係数）を逐次的に決定します。

AR(\(M\))モデルは、現在のデータ \(y_t\) を過去 \(M\) 個のデータの線形結合とノイズで表現します。

\[ y_t = -\sum_{k=1}^{M} a_k y_{t-k} + \epsilon_t \]

ここで、\(a_k\) はAR係数、\(\epsilon_t\) は白色ノイズです。

Burg法では、Levinson-Durbin再帰アルゴリズムを用いて、ARモデルの次数 \(m\) を1から \(\text{mmx}\) まで増やしながらAR係数 \(a_k^{(m)}\) と予測誤差パワー \(P_m\) を更新します。

AR係数 \(a_m^{(m)}\) は、予測誤差を最小化するように以下の式で求められます（プログラム中の pg[m] に相当）。

\[ a_m^{(m)} = - \frac{2 \sum_{i=1}^{N-m} f_{m-1}(i) b_{m-1}(i-1)}{\sum_{i=1}^{N-m} (f_{m-1}(i)^2 + b_{m-1}(i-1)^2)} \]

ここで、\(f_{m-1}(i)\) と \(b_{m-1}(i-1)\) はそれぞれ前向きおよび後ろ向き予測誤差です。プログラム中の sn と sd はこれらの誤差の積和・二乗和に対応します。

予測誤差パワー \(P_m\) は、以下の式で更新されます（プログラム中の pm）。

\[ P_m = P_{m-1} (1 - (a_m^{(m)})^2) \]

AR係数 \(a_k^{(m)}\) は、低次モデルの係数 \(a_k^{(m-1)}\) と \(a_m^{(m)}\) を用いて以下の再帰式で更新されます（プログラム中の pg[k]）。

\[ a_k^{(m)} = a_k^{(m-1)} + a_m^{(m)} a_{m-k}^{(m-1)} \quad (1 \le k < m) \]

最終的に、推定されたAR係数 \(a_k\) と予測誤差パワー \(P_M\) を用いて、パワースペクトル \(P(f)\) は以下の式で計算されます（プログラム中の pspec[i]）。

\[ P(f) = \frac{2 \cdot P_M \cdot \Delta t}{|1 + \sum_{j=1}^{M} a_j e^{-i 2\pi f j \Delta t}|^2} \]

ここで、\(\Delta t\) はサンプリング間隔、\(f\) は周波数です。分母の絶対値の二乗は、プログラム中の sumr * sumr + sumi * sumi に対応します。

プログラムの冒頭で x0, x1, nData, WG, IgnoreZero, ConvFunc, mmx, lmax などのパラメータがハードコードされており、必要に応じて直接編集できます。なお、WG, IgnoreZero, ConvFunc の各変数は、現在のコードでは使用されていません。

必要な非標準ライブラリとインストール方法

このプログラムは、以下の非標準ライブラリを使用しています。

  • NumPy: 数値計算を効率的に行うためのライブラリ。配列操作などで使用されます。

  • Matplotlib: グラフ描画のためのライブラリ。結果の可視化に使用されます。

これらのライブラリは、Pythonのパッケージマネージャ pip を使ってインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。

pip install numpy matplotlib

必要な入力ファイル

mem.py は、外部の入力ファイルを必要としません。 プログラム内部の Func(x) 関数が時系列データを生成します。この関数は、以下の3つのサイン波を合成し、さらに少量のランダムノイズを追加しています。

  • 周波数 1.5 Hz, 振幅 1.0

  • 周波数 3.0 Hz, 振幅 0.3

  • 周波数 10.0 Hz, 振幅 0.5

データ点の数 (nData) や時間範囲 (x0, x1) は、プログラム冒頭のグローバル変数で設定されています。

生成される出力ファイル

mem.py は、実行時にファイルとしてデータを保存しません。

以下の内容が標準出力（コンソール）に出力されます。

  • サンプリング情報 (x0, x1, xstep, nData, W(Gauss))

  • 周波数 (f)、パワースペクトル (spc)、予測誤差パワー (fpe)、および自己相関関数 (auto c.) の計算結果を表形式で出力。出力される列は、計算される次数 mmx や lmax の値によって変化します。

また、matplotlib によって以下のグラフが画面に表示されます。

  • 上段: 元の時系列データ (input)。横軸が時間 t、縦軸がデータ値 y。

  • 下段: 最大エントロピー法で推定されたパワースペクトル (mem)。横軸が周波数 f、縦軸がスペクトル強度 mem。

プログラムは、ユーザーがEnterキーを押すまでグラフを表示したまま待機します。

コマンドラインでの使用例 (Usage)

mem.py は、コマンドライン引数を必要とせず、直接実行することができます。

python mem.py

コマンドラインでの具体的な使用例

以下のコマンドを実行することで、プログラムが生成した時系列データのMEMスペクトル解析が実行され、結果が標準出力に表示され、グラフがポップアップウィンドウで表示されます。

python mem.py

実行後、コンソールには以下のような情報が出力されます（具体的な数値は実行環境や内部乱数により多少異なります）。

x = (0, 1.0, 0.00097751953125)
nData= 1024
W(Gauss)= 0.02

f       spc     fpe     auto c.
0.0     0.0     0.0     -0.1030096969888691
0.97751953125   0.00010998595604117062    0.00010322237937966953  -0.08271701548161555
1.9550390625    0.0016489376189993318     0.0015487779313460875   -0.04609825633830495
2.93255859375   0.0014022699264421422     0.0013165183358988514   0.012586884674753068
... (以下、`lmax`で指定された次数まで、f, spc, fpe, auto c. の4列が出力) ...
37.1457421875   0.0006761664177243987     0.000635105221978252    0.0003024847953689459
38.12326171875  0.0001740924976451624     0.0001633534579172288   -0.0002167198756059639
... (以下、`lmax`まで、f, spc, auto c. の3列が出力) ...
49.853515625    0.000003362077795325852   -0.0000013695277561825707
... (以下、nmaxまで、f, spc の2列が出力) ...
Press ENTER to exit>>

同時に、以下に示すような2つのグラフを含むウィンドウが表示されます。 上段には、生成された時系列データが表示され、その変動を確認できます。 下段には、MEMによって計算されたパワースペクトルが表示され、時系列データに含まれる主要な周波数成分（この例では1.5 Hz, 3.0 Hz, 10.0 Hz付近）がピークとして現れていることを視覚的に確認できます。

（ここにグラフのイメージを挿入することはできませんが、上記の説明が内容を示しています。）

ユーザーがコンソールで ENTER キーを押すと、プログラムとグラフウィンドウが終了します。