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
```

## CWTの実装

まずContinuous wavelet transform(CWT)の式は以下である.　参考: [wiki: Continuous Wavelet Transform(CWT)](https://en.wikipedia.org/wiki/Continuous_wavelet_transform)


${\displaystyle X_{w}(a,b)={\frac {1}{|a|^{1/2}}}\int _{-\infty }^{\infty }x(t){\overline {\psi }}\left({\frac {t-b}{a}}\right)\,dt}$

なんかすごそうに感じるが以下のようにとてもシンプルである.

```python

class MotherWavelet(object):

    def __init__(self, wavelet):
        self.wavelet = wavelet

    def __call__(self, wid, a):
        t = np.arange(wid) - int(wid / 2)
        return self.wavelet(t  / a)


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

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

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

以下でその波形の位置をずらしながら畳み込んでいる(比較しているような感じ).

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

最後にこの時順番にｙに追加しているので, サイズ`i`の波との比較結果はy[i]に格納される.

```python
y.append(s / np.abs(a)**0.5)
```

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

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

# 数式との対応

全体像がつかめたところで数式とちゃんと対応してるか見ていく.

${\displaystyle X_{w}(a,b)={\frac {1}{|a|^{1/2}}}\int _{-\infty }^{\infty }x(t){\overline {\psi }}\left({\frac {t-b}{a}}\right)\,dt}$

---
まず以下で`MotherWavelet`クラスのコール文が呼ばれている.

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

このMotherWaveletクラスのコール文はこの部分.

```python
    def __call__(self, wid, a):
        t = np.arange(wid) - int(wid / 2)
        return self.wavelet(t  / a)
```

そしてこの部分が数式の以下の部分を表している.

${\displaystyle {\overline {\psi }}\left({\frac {t-b}{a}}\right)}$


具体的に言うと, -wid/2からwid/2までの`t`の配列を作って全体を係数`a`で割っている.
更に, ここで生成したwaveは配列tの範囲外では0とみなす. ← 結構重要
${-b}$ が無いように見えるが, これについては次で出てくるので一旦無視してほしい.

---

次に`convolve`があり,作成した${{\overline {\psi }}\left({\frac {t-b}{a}}\right)}$をずらしながら元の波形${x(t)}$と掛け合わせてその相和を求めている.
(相和なのでmodeが'same')

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

そして. この時どれだけずらしたのか（バイアス）が先ほど無視した${b}$に相当する.
また, 明らかに有限の範囲で積を求めているが, tの範囲外は先ほど述べたように0であり, 当然積も0になるので問題ない.
なので, ここまでで以下の部分になる.

${\displaystyle \int _{-\infty }^{\infty }x(t){\overline {\psi }}\left({\frac {t-b}{a}}\right)\,dt}$

---

あとは最後に${\frac {1}{|a|^{1/2}}}$を添えて

```python
y.append(s / np.abs(a)**0.5)
```

完成！ρ(・ω・、)

${\displaystyle X_{w}(a,b)={\frac {1}{|a|^{1/2}}}\int _{-\infty }^{\infty }x(t){\overline {\psi }}\left({\frac {t-b}{a}}\right)\,dt}$


### 最後に

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

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


class MotherWavelet(object):

    def __init__(self, wavelet):
        self.wavelet = wavelet

    def __call__(self, wid, a):
        x = (np.arange(wid) - int(wid / 2)) / a
        return self.wavelet(x)


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 cwt(x, mw, A):
    y = []
    for a in A:
        wave = mw(min(10 * a, len(x)), a)
        s = np.convolve(x, wave, mode='same')
        y.append(s / np.abs(a)**0.5)
    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, MotherWavelet(ricker), widths)
print('x(t)')
plt.plot(sig)
plt.show()
print('Mother Wavelet')
# plt.plot(signal.ricker(200, 10))
plt.plot(MotherWavelet(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を実装した.

# 感想

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