Python: ウェーブレット変換の実装：②CWT編

# はじめに

## 目標

本投稿では, Continuous Wavelet Transform(CWT)の理解のためscpyを使わずにCWT変換を行う関数を作成する.
参考: [wiki: Continuous Wavelet Transform(CWT)](https://en.wikipedia.org/wiki/Continuous_wavelet_transform)

## 参考

まず参考にCWTがどのようなものかscpyのcwtとその変換例を紹介する. 参考: [scipy.signal.cwt](https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.signal.cwt.html)
以下は上記のサンプルコードに対して元波形${x(t)}$をMother Waveretの可視化を追加したものである.
Mother Waveretについての説明はおいおい

```python
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

t = np.linspace(-1, 1, 200, endpoint=False)
sig  = np.cos(2 * np.pi * 7 * t) + np.sin(np.pi * 6 * t)
widths = np.arange(1, 31)
cwtmatr = signal.cwt(sig, signal.ricker, widths)
print('x(t)')
plt.plot(sig)
plt.show()
print('Mother Wavelet')
plt.plot(signal.ricker(200, 10))
plt.show()
print('Xw(a, b)')
plt.imshow(cwtmatr, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
           vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
plt.show()
```

以下サンプルコードの実行例 1つ目が元波形, 2つ目がMother Waveret, 3つ目がCWT変換後である.

![サンプル版.png](https://qiita-image-store.s3.amazonaws.com/0/365961/15914e88-f1e7-c05c-85df-e47863ca040d.png)

# CWTの実装とテスト

## Mther Wavelet
CWTにはMther Waveleが必要になるが, これは以前[ウェーブレット変換の実装：①Ricker Waveletの実装編](https://qiita.com/sai-sui/items/10adf84a4e3761706cd0)で実装した以下を用いる(詳細はリンク先にて)

```python
def ricker(t, s=1.0):
    a = 2 * np.sqrt(3 * s) * np.pi**0.25
    b = 1 - (t / s)**2
    c = np.e**(-1 * t**2 / (2 * s**2))
    return a * b * c


def s_ricker(wid, a):
    x = (np.arange(wid) - int(wid / 2)) / a
    return ricker(x)
```

## CWTの実装
CWTによるウェーブレット変換(｀・ω・´)ｷﾘｯとかカッコつけるとすごそうに感じるが.
今回実装したCWT変換以下のようにとてもシンプルである.

```python
def cwt(x, mw, A):
    y = []
    for a in A:
        wave = mw(min(10 * a, len(x)), a)
        y.append(np.convolve(x, wave, mode='same'))
    return np.array(y)
```

具体的には, Mother Waveletから以下の部分でサイズ感`a`に応じた波を作って

```python
wave = mw(min(10 * a, len(x)), a)
```

以下でその波形の位置をずらしながら畳み込んでいる(比較しているような感じ).
この時順番にｙに追加しているので, サイズ`i`の波との比較結果はy[i]に格納される.

```python
y.append(np.convolve(x, wave, mode='same'))
```

雰囲気が伝わるようにfor文におけるa=5の時の例で簡単な画像を作ってみたので参考にしてほしい.
(注意：あくまで画像はイメージです)

![大まかなイメージ.png](https://qiita-image-store.s3.amazonaws.com/0/365961/fa37506f-c688-5f7a-7843-148e0f931dac.png)


実際に最初のサンプルコードと同じ変換が実現できるか確認してみる.

```python
import matplotlib.pyplot as plt
import numpy as np


def ricker(t, s=1.0):
    a = 2 * np.sqrt(3 * s) * np.pi**0.25
    b = 1 - (t / s)**2
    c = np.e**(-1 * t**2 / (2 * s**2))
    return a * b * c


def s_ricker(wid, a):
    x = (np.arange(wid) - int(wid / 2)) / a
    return ricker(x)


def cwt(x, mw, A):
    y = []
    for a in A:
        wave = mw(min(10 * a, len(x)), a)
        y.append(np.convolve(x, wave, mode='same'))
    return np.array(y)


t = np.linspace(-1, 1, 200, endpoint=False)
sig  = np.cos(2 * np.pi * 7 * t) + np.sin(np.pi * 6 * t)
widths = np.arange(1, 31)
# cwtmatr = signal.cwt(sig, signal.ricker, widths)
cwtmatr = cwt(sig, s_ricker, widths)
print('x(t)')
plt.plot(sig)
plt.show()
print('Mother Wavelet')
# plt.plot(signal.ricker(200, 10))
plt.plot(s_ricker(200, 10))
plt.show()
print('Xw(a, b)')
plt.imshow(cwtmatr, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
           vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
plt.show()
```

![自作版.png](https://qiita-image-store.s3.amazonaws.com/0/365961/25ed5b10-7cfa-8652-b345-a14bd3f3838e.png)

# まとめ

scpyのcwtと（ほぼ）同じCWTを実装した.

# 感想

そこそこ理解できてきた気がしなくもない(´∀｀)