正弦波の振幅と位相をカルマンフィルターで求める

## 概要
pykalmanライブラリのカルマンフィルター を使って、正弦波の振幅と位相を求めてみます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/73020/0f7b0244-7b36-06b1-f302-9cb7166edf84.png)


### 他の手法との比較
正弦波の振幅と位相をScipyのカーブフィットで求める
https://qiita.com/nnn112358/items/35cfee57bd1d8f21147b

正弦波の振幅と位相をSciPyのFFTで求める
https://qiita.com/nnn112358/items/dfded784541045f988d4

正弦波の振幅と位相をLSTMで求める
https://qiita.com/nnn112358/items/1c9b077a7a413ee1dd9f

正弦波の振幅と位相をRNNで求める
https://qiita.com/nnn112358/items/63b53fb192b7f10d509a

### カルマンフィルターの説明

カルマンフィルターは、誤差のある観測値を用いて、ある動的システムの状態を推定あるいは制御するための、無限インパルス応答フィルターの一種です。

具体的には、以下の2つの情報を用いて、システムの状態を推定します。

- システムモデル: システムの動作を記述した数学的なモデル
- 観測値: システムの状態を直接測定した値

システムモデルは、システムの状態変化を記述した微分方程式などで表されます。観測値は、センサーなどで測定した値ですが、一般的には誤差を含んでいます。

カルマンフィルターは、これらの情報に基づいて、システムの現在の状態と、システムの次の状態を予測する推定値を計算します。

## Python開発環境のインストール

Python環境として、UbuntuにMinicondaをインストールする手順します。

### MinicondaのLinuxへのインストール手順

1. **Minicondaのインストーラをダウンロード**:
    まず、Minicondaの公式ウェブサイトからインストーラをダウンロードします。以下のコマンドで最新版のMinicondaインストーラをダウンロードします。

    ```bash
    wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
    ```

2. **インストーラの実行**:
    ダウンロードしたスクリプトを実行してMinicondaをインストールします。インストール先を指定する場合は、オプションでディレクトリを指定できます。

    ```bash
    bash Miniconda3-latest-Linux-x86_64.sh
    ```

    インストールプロセスが開始されると、いくつかのプロンプトが表示されます。これに従って進めてください。

    - 使用許諾契約書が表示されたら、スペースキーを押して全文を表示し、`yes`と入力して同意します。
    - インストールディレクトリの確認が表示されるので、デフォルトのままで良ければそのままエンターキーを押します。カスタムディレクトリにインストールしたい場合は、パスを入力してエンターキーを押します。
    - インストールが完了したら、`conda`コマンドを使えるようにするために、環境変数の設定を行うかどうかのプロンプトが表示されます。これも`yes`と入力して進めます。

3. **環境変数の設定**:
    インストールが完了したら、Minicondaのバイナリディレクトリを環境変数`PATH`に追加します。多くの場合、インストールスクリプトが自動的にこれを行いますが、手動で行う場合は以下のコマンドを実行します。

    ```bash
    export PATH=~/miniconda3/bin:$PATH
    ```

    この設定を永続化するために、`~/.bashrc`（または使用しているシェルの設定ファイル）に追加します。

    ```bash
    echo 'export PATH=~/miniconda3/bin:$PATH' >> ~/.bashrc
    source ~/.bashrc
    ```

4. **インストールの確認**:
    `conda`コマンドを実行して、Minicondaのインストールが正しく行われたことを確認します。

    ```bash
    conda --version
    ```

    これでMinicondaのバージョンが表示されれば、インストールは成功です。

5. **Condaの初期設定とアップデート**:
    インストール後、初期設定を行い、Conda自体とデフォルトパッケージをアップデートします。

    ```bash
    conda init
    conda update conda
    ```

これでMinicondaのインストールは完了です。

### 仮想環境の作成

Minicondaのインストールが完了したら、仮想環境を作成します。

 **仮想環境の作成**:

```bash
 conda create -n new_env python=3.8 tensorflow keras matplotlib scikit-learn pykalman
```


 **仮想環境のアクティベート**:

```bash
conda activate new_env
```

 **仮想環境のディアクティベート**:

 仮想環境を終了するには、以下のコマンドを実行します。

```bash
conda deactivate
```


## サンプルソース

以下に、カルマンフィルターを用いて正弦波の振幅と位相を求める手順を示します。この手順では、カルマンフィルターを用いて信号を予測し、その予測結果から振幅と位相を推定します。

```
$ python kalman_filter.py
Original amplitude: 2.0000
Estimated amplitude: 2.1435
Original phase: 0.7854
Estimated phase: 0.0213
```

```kalman_filter.py
import numpy as np
import matplotlib.pyplot as plt
from pykalman import KalmanFilter

# パラメータ設定
N = 1000  # サンプル数
T = 0.1  # サンプリング間隔
t = np.linspace(0.0, N*T, N)  # 時間軸の生成

# 既知の周期と周波数
known_period = 5.0  # 既知の周期
known_freq = 1.0 / known_period  # 周波数の計算

# 元の正弦波を生成
amplitude_original = 2.0  # 元の正弦波の振幅
phase_original = np.pi / 4  # 元の正弦波の位相
y_original = amplitude_original * np.sin(2.0 * np.pi * known_freq * t + phase_original)  # 元の正弦波データ

# ノイズを加える
np.random.seed(42)  # 再現性のためにランダムシードを設定
noise = np.random.normal(0, 0.5, N)  # 正規分布に従うノイズを生成
y_noisy = y_original + noise  # ノイズを加えたデータ

# カルマンフィルターの設定
transition_matrix = [[1, -T * 2 * np.pi * known_freq], [T * 2 * np.pi * known_freq, 1]]
observation_matrix = [[1, 0]]

kf = KalmanFilter(transition_matrices=transition_matrix,
                  observation_matrices=observation_matrix,
                  initial_state_mean=[0, 1],
                  initial_state_covariance=np.eye(2),
                  observation_covariance=1,
                  transition_covariance=np.eye(2) * 0.1)

# カルマンフィルターを適用して状態を推定
states_pred = kf.em(y_noisy).filter(y_noisy)[0]

# 振幅と位相の推定
estimated_amplitude = np.sqrt(states_pred[:, 0]**2 + states_pred[:, 1]**2)
estimated_phase = np.arctan2(states_pred[:, 1], states_pred[:, 0])

# 平均振幅と位相を求める
mean_amplitude = np.mean(estimated_amplitude)
mean_phase = np.mean(estimated_phase)

# 結果の表示
print(f"Original amplitude: {amplitude_original:.4f}")
print(f"Estimated amplitude: {mean_amplitude:.4f}")
print(f"Original phase: {phase_original:.4f}")
print(f"Estimated phase: {mean_phase:.4f}")

# フィッティングした正弦波を生成
y_fit = mean_amplitude * np.sin(2.0 * np.pi * known_freq * t + mean_phase)

# グラフ描画
plt.figure(figsize=(14, 8))

plt.plot(t, y_noisy, label='Noisy Data', alpha=0.5)
plt.plot(t, y_original, label='Original Sine Wave', linewidth=2)
plt.plot(t, y_fit, label='Fitted Sine Wave (Kalman)', linestyle='--')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()

plt.tight_layout()
plt.show()

```