Pythonで連続ウェーブレット変換(scipy, mlpy, swan)

はじめに

そもそもウェーブレット変換ってなに？　という方はこの記事を読まないと思うのでそのあたりの説明ははしょります．が，いざpythonで連続ウェーブレット変換（Continuous WaveletTransform：cwt）しようとした際どハマりしたので備忘録を残しておきます．

pythonでcwtあれこれ

今回ウェーブレット変換をかける時系列データとして，単純で結果が見やすいサンプルデータを生成します．

import matplotlib.pyplot as plt
import numpy as np

x = np.arange(0, 20, 0.01)
y = np.sin(2 * np.pi * 1 * x) * 2 + np.sin(2 * np.pi * 2 * x) * 2  + np.cos(2 * np.pi * 10 * x)   

plt.plot(x, y)
plt.show()

上図のように3つの波形を重ね合わせたものを用います．

scipy.signal.cwt()

scipyの公式サイトにサンプルが記載されているのでそれを使用してみます．

from scipy import signal
widths = np.arange(1, 31)
cwtmatr = signal.cwt(y, signal.ricker, widths)
plt.imshow(cwtmatr, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
            vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())  

例では，マザーウェーブレットにsignal.ricker，メキシカンハットを利用しています．
ただしこの結果，横軸は時間でいいとしても縦軸が何かよくわかりません．
StackOverflowにも質問があがっていますが，どうやらスケールではないみたい... なのでマザーウェーブレットにモルレーを用いてスケールから擬似周波数に変換するなんてことも難しそう...

PyWavelets

pythonでwavelet変換を行うライブラリとして，PyWaveletsというものが公開されています．Qittaにもいくつかこちらのライブラリを使った記事が上がっていますので詳しい使い方はそちらをご参照いただければと思います．

[参考]
HirofumiYashimaさま
【 Python で連続ウェーブレット変換 】 PyWavelets モジュール の「公式マニュアル」に書いていない 「マザーウェーブレット」略語 一覧情報 記載場（https://qiita.com/HirofumiYashima/items/3e89fe1e79f8079b41b6）

さて，こちらを利用してcwtを試そうとしたのですが，こちらのpywt.cwt()関数も，出力されているものとスケールの対応がどうなっているのかわかりませんでした．

import pywt
widths = np.arange(1, 31)
cwtmatr, freqs = pywt.cwt(y, widths, 'mexh')
plt.imshow(cwtmatr, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
            vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())  

検索して出てくるものは主に上記2つであり，八方塞がりだったのですが，調べているうちに機械学習用ライブラリであるmlpyにもcwt用関数が入っているということがわかったのでそちらを試してみることにしました．

mlpy

ただ，このmlpyは如何せん古いライブラリで，更新も2012年で止まってしまっています．インストールに手間取ったので備忘録を残しておきます．

mlpyのインストール

環境

• macOS Sierra (10.12.4)
• Python 3.5.2

インストール

pipではうまく入れられず色々問題が多そうだったので，ソースからインストールします．
まず，リンクから最新版のソースを落としてきます(mlpy-3.5.0.tar.gz)．
解凍したのち，mlpy-3.5.0に移動してpython setup.py installとやると色々怒られます．
glsが必要とのことなのでbrew install glsでインストールします(参考)．

また，ソースコードを何箇所か修正する必要があるようなので，参考サイトをみつつ修正します．

[参考サイト]
soyiharuさま
mlpyをpython3.3(64bit)でビルド(windows 64bit)

これに加え，連続ウェーブレット変換をかける場合は，エラー回避のため以下を修正します．

修正前

mlpy/wavelet/continuous.py（161~165行目）

    J = floor(dj**-1 * log2((N * dt) / s0))

    s = empty(J + 1)

    for i in range(s.shape[0]):
        s[i] = s0 * 2**(i * dj)

修正後

mlpy/wavelet/continuous.py（161~165行目）

    J = floor(dj**-1 * log2((N * dt) / s0))

    s = empty(int(J) + 1)

    for i in range(s.shape[0]):
        s[i] = s0 * 2**(i * dj)

ここまで終わったら，

 python setup.py build_ext --include-dirs=/path/to/gsl/include --rpath=/path/to/gsl/lib install                  

でインストールします．brewでインストールした場合，includeファイルまではパスを通してくれないようなので，上記のような指定をする必要があります．

###mlpyで連続ウェーブレット変換

ここまですんだら後はコードを書いていくだけです．作成にあたりこちらを参考にしました．
こちらのリンクにあるように，obspyにもcwt機能があるようですが，データを一度mseed形式に変換する必要がありそうだったので，今回は使用を見送りました．

import mlpy

omega0 = 8
wavelet_fct = "morlet"
# (1)スケール設定
scales = mlpy.wavelet.autoscales(N=len(x), dt=0.01, dj=0.05,
                                 wf=wavelet_fct, p=omega0)
# (2)連続ウェーブレット変換
spec = mlpy.wavelet.cwt(y, dt=0.01, scales=scales,
                        wf=wavelet_fct, p=omega0)
# (3)スケールから擬似周波数を算出
freq = (omega0 + np.sqrt(2.0 + omega0 ** 2)) / (4 * np.pi * scales[1:])

plt.rcParams['figure.figsize'] = (20, 6)
fig = plt.figure()
ax1 = fig.add_axes([0.1, 0.75, 0.7, 0.2])
ax2 = fig.add_axes([0.1, 0.1, 0.7, 0.60], sharex=ax1)
ax3 = fig.add_axes([0.83, 0.1, 0.03, 0.6])

ax1.plot(x, y, 'k')

# (4)振幅を取るためにspecの絶対値をとり，縦軸をfreqに
img = ax2.imshow(np.abs(spec), extent=[0, 20, freq[-1], freq[0]],
                 aspect='auto', interpolation='nearest')
# Hackish way to overlay a logarithmic scale over a linearly scaled image.
twin_ax = ax2.twinx()
twin_ax.set_yscale('log')
twin_ax.set_xlim(0, 20)
twin_ax.set_ylim(freq[-1], freq[0])
ax2.tick_params(which='both', labelleft=False, left=False)
twin_ax.tick_params(which='both', labelleft=True, left=True, labelright=False)

fig.colorbar(img, cax=ax3)

plt.show()

関数の詳しい使い方はmlpyの公式サイトを参照いただければ良いかと思います．

(1)はmlpy.wavelet.autoscalesでスケールの設定を行なっており，引数として

N	サンプル数
dt	データのタイムステップ(例では0.01s)
dj	スケール解像度(小さい方が高解像度)
wf	マザーウェーブレット(morlet, paul, dog)
p	マザーウェーブレットのパラメータ(morletはω0, paul, dogならorder)

が指定されます．返り値として，スケールのnp.arrayが得られます．

(2)はmlpy.wavelet.cwtで連続ウェーブレット変換を行なっており，引数として

x	データ(一次元のnp.array)
dt	データのタイムステップ(autoscaleと揃える)
scales	autoscalesで得られたスケールの配列
wf	マザーウェーブレット(autoscaleと揃える)
p	マザーウェーブレットのパラメータ(autoscaleと揃える)

が指定されます．

(3)では，以下の式に基づいてスケールから擬似周波数を算出しています．

f' = \frac{\omega_0 + \sqrt{2 + \omega_0^2}}{4\pi s} 

以上の結果，

上図のように，縦軸で，1,2,10あたりに高い数値のバンドが現れているので，縦軸が正しくFreqencyとして出力されていることがわかります．

swan

swanというwavelet変換用ライブラリは最新更新が　２０１７年2月と比較的新しいのでこちらも試してみます．

こちらのインストールはpipで問題なく入りました．

from swan import pycwt

Fs = 1/0.01
omega0 = 8
# (1)　Freqを指定してcwt
freqs=np.arange(0.1,20,0.025)
r=pycwt.cwt_f(y,freqs,Fs,pycwt.Morlet(omega0))
rr=np.abs(r)

plt.rcParams['figure.figsize'] = (20, 6)
fig = plt.figure()
ax1 = fig.add_axes([0.1, 0.75, 0.7, 0.2])
ax2 = fig.add_axes([0.1, 0.1, 0.7, 0.60], sharex=ax1)
ax3 = fig.add_axes([0.83, 0.1, 0.03, 0.6])

ax1.plot(x, y, 'k')

img = ax2.imshow(np.flipud(rr), extent=[0, 20, freqs[0], freqs[-1]],
                 aspect='auto', interpolation='nearest')

fig.colorbar(img, cax=ax3)

plt.show()

swanでは(1)にあるように， Freqを自分で指定しcwtにかけることができるようです．
出力は以下のようになり，

縦軸はlogスケールになっていませんが，正しいバンドが出力されていることがみてとれます．

おわりに

pythonでの連続ウェーブレット変換が思った以上に厄介だったので簡単にまとめました．
FFTに比べるとなかなか情報が少なく，Rやmatlabを使用する例が多くみられるので，pythonでの処理はあまり一般的ではないのかもしれませんが，参考になればと思います．
scipy.signalやPyWaveletsにもmlpyのautoscaleのような機能があるのかもしれませんが，こちらの調査不足で見つけることが叶わなかったので，もしご存知の方は教えいただけると幸いです．

参考

• 連続ウェーブレット変換に使うマザーウェーブレット色々: Morlet, Paul, DOG（http://r9y9.github.io/blog/2014/06/01/continuouos-wavelet-transform-types/）
• Torrence, C. and G.P. Compo “A Practical Guide to Wavelet Analysis”, Bull. Am. Meteorol. Soc., 79, 61–78, 1998.（http://paos.colorado.edu/research/wavelets/bams_79_01_0061.pdf）
• ウェーブレット変換の基礎と応用事例：連続ウェーブレット変換を中心に（https://www.slideshare.net/ryosuketachibana12/ss-42388444）