はじめに
目標
本投稿では, 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
CWTの実装
まずContinuous wavelet transform(CWT)の式は以下である. 参考: wiki: Continuous Wavelet Transform(CWT)
${\displaystyle X_{w}(a,b)={\frac {1}{|a|^{1/2}}}\int _{-\infty }^{\infty }x(t){\overline {\psi }}\left({\frac {t-b}{a}}\right),dt}$
なんかすごそうに感じるが以下のようにとてもシンプルである.
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に応じた波を作って
wave = mw(min(10 * a, len(x)), a)
以下でその波形の位置をずらしながら畳み込んでいる(比較しているような感じ).
s = np.convolve(x, wave, mode='same')
最後にこの時順番にｙに追加しているので, サイズiの波との比較結果はy[i]に格納される.
y.append(s / np.abs(a)**0.5)
雰囲気が伝わるようにfor文におけるa=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}$
まず以下でMotherWaveletクラスのコール文が呼ばれている.
wave = mw(min(10 * a, len(x)), a)
このMotherWaveletクラスのコール文はこの部分.
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')
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}}}$を添えて
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}$
最後に
実際に最初のサンプルコードと同じ変換が実現できるか確認してみる.
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()
まとめ
scpyのcwtと（ほぼ）同じCWTを実装した.
感想
結構理解できてきた気がしなくもない(´∀｀)