概要

SciPyライブラリのFFT（高速フーリエ変換）を使って、正弦波の振幅と位相を求めてみます。

他の手法との比較

正弦波の振幅と位相を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

FFTの説明

scipy.fftpack.fft(x, n=None, axis=-1, overwrite_x=False)
パラメータ
    x: 入力データ。時間領域の信号データ。1次元または多次元の配列。
    n: 変換するポイントの数。デフォルトでは、入力データの長さに設定されます。指定した場合、入力データがゼロパディングまたは切り捨てられます。
    axis: 変換する軸。デフォルトは最後の軸。
    overwrite_x: 入力データを上書きして計算を行うかどうか。メモリ使用量を削減するために使用されますが、元のデータが変更されます。
戻り値
    out: 周波数領域の複素数配列。入力データと同じ形状。

fft関数は、離散フーリエ変換 (Discrete Fourier Transform, DFT) を計算するための効率的なアルゴリズムである高速フーリエ変換 (Fast Fourier Transform, FFT) を実装したものです。SciPyライブラリのscipy.fftpackモジュールで提供されています。この関数で、時間領域の信号を周波数領域に変換し、信号の周波数成分を解析することができます。

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

Python環境として、Minicondaをインストールする手順を説明します。
WSLにUbuntu22.04 をインストールして使用しました。

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

Minicondaのインストーラをダウンロード:
まず、Minicondaの公式ウェブサイトからインストーラをダウンロードします。以下のコマンドで最新版のMinicondaインストーラをダウンロードします。
```
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
```
インストーラの実行:
ダウンロードしたスクリプトを実行してMinicondaをインストールします。インストール先を指定する場合は、オプションでディレクトリを指定できます。
```
bash Miniconda3-latest-Linux-x86_64.sh
```
インストールプロセスが開始されると、いくつかのプロンプトが表示されます。これに従って進めてください。
- 使用許諾契約書が表示されたら、スペースキーを押して全文を表示し、yesと入力して同意します。
- インストールディレクトリの確認が表示されるので、デフォルトのままで良ければそのままエンターキーを押します。カスタムディレクトリにインストールしたい場合は、パスを入力してエンターキーを押します。
- インストールが完了したら、condaコマンドを使えるようにするために、環境変数の設定を行うかどうかのプロンプトが表示されます。これもyesと入力して進めます。
環境変数の設定:
インストールが完了したら、Minicondaのバイナリディレクトリを環境変数PATHに追加します。多くの場合、インストールスクリプトが自動的にこれを行いますが、手動で行う場合は以下のコマンドを実行します。
```
export PATH=~/miniconda3/bin:$PATH
```
この設定を永続化するために、~/.bashrc（または使用しているシェルの設定ファイル）に追加します。
```
echo 'export PATH=~/miniconda3/bin:$PATH' >> ~/.bashrc
source ~/.bashrc
```
インストールの確認:
condaコマンドを実行して、Minicondaのインストールが正しく行われたことを確認します。
```
conda --version
```
これでMinicondaのバージョンが表示されれば、インストールは成功です。
Condaの初期設定とアップデート:
インストール後、初期設定を行い、Conda自体とデフォルトパッケージをアップデートします。
```
conda init
conda update conda
```

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

仮想環境の作成

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

仮想環境の作成:

 conda create -n new_env python=3.8 matplotlib scikit-learn

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

conda activate new_env

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

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

conda deactivate

サンプルソース

正弦波データにノイズを加え、そのデータに対してFFT関数を使用して最適な正弦波のパラメータを推定し、結果をグラフに表示してみます。

$ python fft.py
Original amplitude: 2.0000
Fitted amplitude: 2.0313
Original phase: 0.7854
Fitted phase: 0.7442 radians
Fitted frequency: 0.2000 Hz

fft.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, fftfreq

# パラメータ設定
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  # ノイズを加えたデータ

# FFTを実行
y_fft = fft(y_noisy)
freqs = fftfreq(N, T)

# 正弦波の振幅と位相を求める
idx = np.argmax(np.abs(y_fft[:N//2]))  # 最大ピークを持つ周波数のインデックスを見つける
fft_freq = freqs[idx]
fft_amplitude = 2.0 / N * np.abs(y_fft[idx])  # 振幅を求める
fft_phase = -np.angle(y_fft[idx])  # 位相を求める

# 結果表示
print(f"Original amplitude: {amplitude_original:.4f}")
print(f"Fitted amplitude: {fft_amplitude:.4f}")
print(f"Original phase: {phase_original:.4f}")
print(f"Fitted phase: {fft_phase:.4f} radians")
print(f"Fitted frequency: {fft_freq:.4f} Hz")

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

# グラフ描画
plt.figure(figsize=(12, 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', linestyle='--')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Sine Wave Fitting with FFT')
plt.legend()

plt.tight_layout()
plt.show()