8-5. CWT_Acc.py – BiP Help

(time, acceleration_longitudinal,acceleration_lateral,acceleration_dorso_ventral) を含むデータ（data.csv）*読み込み、３軸の加速度を元にODBAを計算し、それを元に以下のような連続ウェーブレット変換による周波数スペクトルを表示するツールです。

*BiPでOpenデータとなっているKatsufumi Sato (Atmosphere and Ocean Research Institute, University of Tokyo) 提供のオオミズナギドリのデータ (title: 9B41870_TS-AxyTrek_Movebank_YNo.6_release20210824)の一部を使用しています。

下記のコードをテキストエディタにコピーし、CWT_Acc.pyという名前で保存します。
scales = np.arange(1, 256)：スケール範囲を0-8Hzに設定しています。
sampling_period=1/20：data.csvのサンプリング周波数（20Hz）に設定しているので、適宜調整してください。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pywt
import matplotlib.dates as mdates

# CSVファイルの読み込み
data = pd.read_csv('data.csv')

# 日時情報をpandasのdatetime型に変換
data['time'] = pd.to_datetime(data['time'])

# ODBAの計算
data['ODBA'] = (
    data['acceleration_longitudinal'].abs() +
    data['acceleration_lateral'].abs() +
    data['acceleration_dorso_ventral'].abs() - 9.8
)

# ウェーブレット変換の設定
odba = data['ODBA'].values
time = data['time']

# 使用するウェーブレットとスケールの設定（具体的なパラメータを指定）
wavelet = 'cmor1.5-6.0'  # 複素モルレウェーブレット（帯域幅1.5、中心周波数6.0）
scales = np.arange(1, 256)  # スケール範囲（調整して周波数を0-8Hzに合わせる）

# 連続ウェーブレット変換
coefficients, frequencies = pywt.cwt(odba, scales, wavelet, sampling_period=1/20)  # サンプリング周波数20Hzと仮定

# グラフの作成
fig = plt.figure(figsize=(10, 8))

# ODBAの生波形のプロット（上部30%）
ax1 = fig.add_axes([0.1, 0.7, 0.8, 0.2])  # [左, 下, 幅, 高さ]
ax1.plot(time, odba, color='red', label='ODBA')
ax1.set_ylabel('ODBA (m/s²)')
ax1.set_title('ODBA and Continuous Wavelet Transform Spectrum')
ax1.grid(True)
ax1.legend(loc='upper right')

# ウェーブレットスペクトルのプロット（下部70%）
ax2 = fig.add_axes([0.1, 0.1, 0.8, 0.55])  # [左, 下, 幅, 高さ]
im = ax2.imshow(
    np.abs(coefficients), 
    extent=[mdates.date2num(time.min()), mdates.date2num(time.max()), frequencies.min(), frequencies.max()],
    cmap='jet', 
    aspect='auto',
    origin='lower'
)

# 周波数の表示範囲を設定（0-8 Hz）
ax2.set_ylim(0, 8)

# 時間を適切なフォーマットで表示
ax2.xaxis_date()
ax2.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))

# 軸ラベルの設定
ax2.set_xlabel('Time')
ax2.set_ylabel('Frequency (Hz)')

# カラーバーの追加
fig.colorbar(im, ax=ax2, label='Amplitude')

# グラフの表示
plt.tight_layout()
plt.show()

pip install matplotlib numpy pywavelets

実行方法

使用するPython環境に必要なライブラリ（pandas, matplotlib）をインストールしてください。

pip install matplotlib numpy pywavelets

コマンドラインでスクリプトを保存したディレクトリに移動し、以下のコマンドを実行します。

python CWT_Acc.py