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

# はじめに

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

# pythonでcwtあれこれ

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

```py3
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()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/59864/af314af7-1ffc-ffeb-fd96-d74e1445f1cc.png)

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

## scipy.signal.cwt()

[scipyの公式サイト](https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.signal.cwt.html)にサンプルが記載されているのでそれを使用してみます．

```py3
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())  
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/59864/7cf113a0-8819-d376-6b6d-2a82960e250d.png)

例では，マザーウェーブレットにsignal.ricker，メキシカンハットを利用しています．
ただしこの結果，横軸は時間でいいとしても縦軸が何かよくわかりません．
[StackOverflowにも質問があがっています](https://stackoverflow.com/questions/28822327/units-of-widths-argument-to-scipy-signal-cwt-function)が，どうやらスケールではないみたい... なのでマザーウェーブレットにモルレーを用いてスケールから擬似周波数に変換するなんてことも難しそう...

## PyWavelets

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

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

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

```py3
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())  
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/59864/7cf113a0-8819-d376-6b6d-2a82960e250d.png)

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

## mlpy

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

### mlpyのインストール

#### 環境

- macOS Sierra (10.12.4)
- Python 3.5.2

#### インストール

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

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

[参考サイト]
soyiharuさま
[mlpyをpython3.3(64bit)でビルド(windows 64bit)](https://qiita.com/soyiharu/items/8db8a50bd04936686c3a)

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

修正前

```py3: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)
```

修正後

```py3: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で連続ウェーブレット変換

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

```py3
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の公式サイト](http://mlpy.sourceforge.net/docs/3.2/wavelet.html)を参照いただければ良いかと思います．

(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)では，以下の式に基づいてスケールから擬似周波数を算出しています．

```math

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

```

以上の結果，

![image.png](https://qiita-image-store.s3.amazonaws.com/0/59864/fc25d951-157f-4b4a-49f3-20fea7e08c34.png)

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

## swan

[swan](https://pypi.python.org/pypi/swan)というwavelet変換用ライブラリは最新更新が　２０１７年2月と比較的新しいのでこちらも試してみます．

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

```py3
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にかけることができるようです．
出力は以下のようになり，

![image.png](https://qiita-image-store.s3.amazonaws.com/0/59864/ec3ed6e8-e960-11e7-fcc1-920d88b3af5e.png)

縦軸は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）