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

はじめに

目標

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

参考

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

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変換後である.

CWTの実装とテスト

Mther Wavelet

CWTにはMther Waveleが必要になるが, これは以前ウェーブレット変換の実装：①Ricker Waveletの実装編で実装した以下を用いる(詳細はリンク先にて)

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変換以下のようにとてもシンプルである.

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に応じた波を作って

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

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

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

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

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

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

まとめ

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

感想

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