工学と数学の間を漂う~

数学が好きな工学徒がメモ代わりに知識を垂れ流します

【Python】Pywaveletsモジュールで連続ウェーブレット変換をしよう

前置き

 連続ウェーブレット変換(以下CWT)を試す必要が出てきたのですが、時間周波数解析の手法としては短時間フーリエ変換(窓付きフーリエ変換)がメジャーであり、ウェーブレット変換の歴史が浅いことが原因で情報が少なく苦戦しました。所説ありますが、短時間フーリエ変換は窓関数の設定が必要でそれらが面倒だということで、窓の調整を自動的に行ってくれるウェーブレット変換が生まれたと言われており、そんな性質を持つウェーブレット変換なので本来、波形データをそのままウェーブレット変換にかければお手軽に時間周波数解析が行えるのがメリットなのですが、情報が少ないのでお手軽な時間周波数解析とは程遠い状況でした。 そこで、これを読めば詳細は置いといてとりあえずPythonでCWTができるというような資料を残そうと思いこの記事を執筆するに至りました。
主旨が「とりあえずCWTをしてみよう」なので実装から天下り的に解説していく構成にしてみました。

Pywaveletsを使ったCWT

 さっそくコード書いて行きます。PythonのPywaveletsモジュールを使ってCWT行いますので、Python?モジュール?なにそれって人は環境構築からモジュールの導入までかるーく調べて準備をしてから以下のコードを動かしてください。 振幅と周波数が時間とともに大きくなるテスト波形を生成し、そのCWTを計算するコードです。

import pywt
import numpy as np
import matplotlib.pyplot as plt

time = np.arange(0, 4, 0.01)
signal = time*np.sin(time**2*np.pi)

dt = time[1]-time[0]
n = len(time)
freq_fft = np.delete(np.fft.rfftfreq(n, dt), 0)

wavelet = "mexh"
scale = pywt.frequency2scale(wavelet, freq_fft*dt)
cwt, freq_cwt = pywt.cwt(signal, wavelet=wavelet, scales=scale)

fig = plt.figure(figsize=(8, 10))
ax1 = fig.add_subplot(211)
ax1.set_xlim(0, 4)
ax1.set_xlabel("Time [s]")
ax1.set_ylabel("Value")
ax1.plot(time, signal)

ax2 = fig.add_subplot(212)
ax2.set_xlabel("Time [s]")
ax2.set_ylim(min(freq_cwt/dt), 30)
ax2.set_ylabel("Frequencies [Hz]")
ax2.pcolormesh(time, freq_cwt/dt, cwt)
plt.show()

 実行結果は以下の通りです。時間とともに周波数が大きくなり、振幅の増大に合わせて色が濃くなっていることが観察できるため時間周波数解析に成功していそうです。

Pywaveletsを使う理由について

 とりあえずCWTはできましたが、わざわざPywaveletsを使用するのはなぜなのでしょうか。それを知るためにはCWTの計算手順を知る必要があります。大雑把な手順は以下の通りです。

  1. 波形データに含まれる周波数を把握する
  2. 波形データに含まれる周波数からスケールを決定
  3. 波形データとスケールからCWTを計算

Pywaveletsでは手順の2に該当する計算をpywt.frequency2scaleにより計算し、3に該当する計算をpywt.cwtで計算しています。また、ウェーブレット変換ではマザーウェーブレット関数を選択する必要があり、短時間フーリエ変換でいうところの窓関数のようにどの関数を使うかによって計算結果が変わってきます。このモジュールではマザーウェーブレット関数がいくつか用意されており、文字列で指定するだけで選択できます。したがって、波形データと波形データの持つ周波数さえわかれば簡単にCWTが計算できます。さらに、周波数はFFTアルゴリズムで簡単に計算ができますので、お手軽という訳です。
纏めると、Pywaveletsはマザーウェーブレット関数の設定、スケールの決定、CWTの計算を手軽に行えるので、主旨にぴったりということです。

そもそもCWTって何?

 Pywaveletsの仕様の説明の中で、スケールだったり、マザーウェーブレット関数だったり初めて聞く用語がいくつかあったと思います。 それらを説明するためには数学的な話をする必要がありますので、ここで説明していきます。(厳密な議論は出来ないのでイメージだけ掴んでいただけると幸いです。)
 まず、フーリエ逆変換は異なる周波数を持った三角関数を無限に足すことで時間関数を再現する手法だというのはご存知だと思います。周波数は連続的な値を取って足されていきますので、以下のように積分により定義されます。

積分なので通常級数のように考えることは無いのですが、三角関数の周波数が連続的に変化しながら足されている級数をイメージしたとき、各項の係数(=振幅)を与えるのがフーリエ変換により計算される関数で、以下のように定義されています。

ここで重要なのはフーリエ逆変換は異なる周波数を持った三角関数の足し合わせで時間関数を再現する手法で、各項の係数であるフーリエ変換は周波数と振幅の関係を与える関数だということです。つまり、ある周波数に対して返される振幅が大きな値を取るとき時間関数においてその周波数をもつ三角関数が多く含まれているという情報を与えてくれます。これがいわゆるフーリエ変換による周波数解析です。

 フーリエ変換と同じような論理展開でウェーブレット変換について説明していきます。ウェーブレット逆変換は三角関数ではなく、マザーウェーブレット関数と呼ばれる波の一部を切り取ったような概形を持つ関数に拡大縮小・並行移動の操作を施したウェーブレット関数を無限に足し合わせることで時間関数を再現します。マザーウェーブレット関数とウェーブレット関数のイメージが難しいと思いますので、ものすごく単純な例を作成してみました。

この例では左辺のような概形を持つ時間関数を異なる拡大縮小と平行移動を施した2つのウェーブレット関数の和で表しています。このようなウェーブレット関数の足し合わせで時間関数が再現できそうなイメージを持てましたでしょうか。拡大縮小を表すパラメータと平行移動を表すパラメータの2つが存在しているので、ウェーブレット関数の足し合わせで時間関数を再現する手法は以下の重積分により定義されます。

C_ψはマザーウェーブレット関数によって決まる定数で、aが拡大縮小を表すパラメータ(=スケール)、bが平行移動を表すパラメータ(=シフト)です。また、W(a, b)という関数が重積分級数でイメージしたときのウェーブレット関数の各項の係数にあたり、ウェーブレット変換と呼ばれ以下の式で定義されています。

フーリエ変換と同様にウェーブレット変換により得られる関数はスケールとシフトと各項の係数の関係を与えます。しかし、このままでは周波数とも時間とも結びつかず、時間周波数解析になりません。スケールとシフトについて以下の換算・解釈を与えることで真価を発揮するようになります。

  • スケールは時間関数の持つ周波数と反比例関係にあり、後述する関係式によって相互に換算が可能

  • シフトは時間軸上の並行移動を表しているので、シフトの値を時間と解釈する

これらを用いることで、W(a, b)は周波数と時間と各項の係数を結びつける関数となるのです。また、フーリエ変換と同じく、あるa(=周波数), b(=時間)を持つウェーブレット関数の係数が大きな値を取るほど、時間関数においてある時間にある周波数をもつウェーブレット関数が多く含まれているという情報を与えてくれますので、時間周波数解析が成立するわけです。

ウェーブレット変換をPCで処理する

 数学的な定義やイメージについて述べましたが、実際のところ波形データはオシロスコープやマイク等により取得されるため時間関数が離散的になります。したがって、前述のウェーブレット変換の定義式を使ってそのまま計算することはできず、積分の離散化が必要になります。
 話は変わりますが、フーリエ変換フーリエ逆変換も積分によって定義されているものの、離散的な波形データ対して計算できますよね。それは積分を離散化した上で高速計算が可能なFFT, IFFTアルゴリズムが存在しているためです。ということはウェーブレット変換をうまく変形してフーリエ変換を使った形に書き換えてやることができれば、積分の離散化に触れることなく計算が行えるということです。以下の式を見てください。 実はウェーブレット変換は畳み込み積分になっているため、畳み込みの定理を使って三行目のようにフーリエ変換フーリエ逆変換を使った形に書き換えることができます。したがって、FFTアルゴリズムとIFFTアルゴリズムを使ってウェーブレット変換をしてやればいいという訳です。しかし、まだ問題があります。それはFFTというかコンピュータは文字を含んだまま計算ができないということです。三行目を見るとψの中に-aが入ってしまっています。これを解決するためには、aに具体的な値を入れ計算三行目の計算を何度も繰り返す必要があります。ですが、aに具体的な値を入れながら計算を行うためにはaにどんな値を取らせるか決めて置く必要があるのと同義です。そこで活躍するのが以下の周波数とスケールの換算式です。

予め波形データに対しフーリエ変換を行い、含まれる周波数を把握しておき、上記の式を使って周波数をスケールaに換算することでaの取り方を定めることができます。
※この周波数とスケールの換算式はマザーウェーブレット関数にメキシカンハット関数を使用した場合にのみ有効です。

改めてコードを眺める

 これまでの説明を踏まえてはじめに載せたコードを見ていきます。

dt = time[1]-time[0]
n = len(time)
freq_fft = np.delete(np.fft.rfftfreq(n, dt), 0)

まず、こちらは波形データのサンプリング間隔とデータ点数を求め、FFTにより周波数を求めています。
※np.deleteは配列の1番目の要素である0 Hzを取り除くためです。換算式が反比例なので、0 Hzを取り除かないとエラーになります。

wavelet = "mexh"
scale = pywt.frequency2scale(wavelet, freq_fft*dt)
cwt, freq_cwt = pywt.cwt(signal, wavelet=wavelet, scales=scale)

こちらは、マザーウェーブレット関数にメキシカンハット関数を指定し、pywt.frequency2scaleを使って周波数からスケールの取り方を決定しています。
その後、pywt.cwtを使って波形データと周波数からCWTを計算しています。
※周波数にサンプリング間隔dtを乗じていますが、Pywaveletsの仕様によるものです。プロット時にはdtで除して周波数を取り出します。

あとがき

 Pywaveletsを使ったCWTの説明を取り扱ったサイトは非常に少ないため、当初は自力でのCWTの実装を目指していました。しかし、Pythonを勉強しながらの実装だったので、挫折し再びPywaveletsを使ったCWTに挑戦しました。自力での実装を目指している間に最低限知っておくべきウェーブレット変換の基礎知識が身についたこともあり、改めてリファレンスを読むとモジュールに用意されている関数の使い道がすんなりと入ってきてなんとかCWTが出来たというのが正直な背景です。ですが、本来ウェーブレット変換は手軽に出来ることを目標に生み出されたものなので、情報の少なさと知名度の低さ故に活躍の機会を失っているようにも感じました。ですので、とりあえず最初のコードをコピーしてtimeとsignalに任意のデータを入れてやることでCWTができるような記事を書こうと思い至った訳です。数学の知識もPythonの知識もまだまだ未熟で、説明も見苦しいものかと思いますが、誰かの御役に立てましたら幸いです。また、質問やご指摘等有りましたらコメントしていただけるとありがたいです。

おまけ

 オシロスコープで波形を取得した場合にtimeとsignalに時間情報と波形情報を割り当てるコード
csvから読み取ることを想定していて、time_columnsに時間情報の列番号を, signal_columnsに波形情報の列番号を指定することでtimeとsignalに 時間情報と波形情報を割り当てることができます。
また、divisionの値はデータ点数を1/数値に削減するための変数です。オシロスコープで取得した波形データのデータ点数が多く、メモリエラーになるときにお試しください。不要であれば削除してつかってください。

import pandas as pd

path = "ここに波形データのパスを指定"
delimiter = "\t"
time_columns = 0
signal_columns = 1
division = 128

data = pd.read_csv(path, delimiter=delimiter)
sampdata = data.iloc[0:0]
for i in np.arange(0, len(data), division):
    sampdata = pd.concat([sampdata, data[i:i+1]])
    
time = sampdata.to_numpy().T[time_columns]
signal = sampdata.to_numpy().T[signal_columns]