はじめに
目標
本投稿では, 以下を参考にContinuous Wavelet Transform(CWT)の逆変換を実装する.
参考: wiki: Continuous Wavelet Transform(CWT)
実装
とりあえずwikiより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}$
${\displaystyle x(t)=C_{\psi }^{{-1}}\int_{{-\infty }}^{{\infty}}\int_{{-\infty }}^{{\infty}}X_{w}(a,b)\frac{1}{|a|^{{1/2}}}{\eta}\left({\frac {t-b}{a}}\right),db\ {\frac {da}{a^{2}}}}$
数式がわかれば後は実装するだけ！できるだけ数式まんまになるようにcwtとicwtを作ってみました.
db, daは中に入れた方が数式まんまになるのですが, 演算量が減るので外に出してます
※変換にはψが必要で, 今回は以前つくったrickerを使う.
import numpy as np
import matplotlib.pyplot as plt
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, p, t, a):
dt = t[1] - t[0]
A = a.reshape(-1, 1, 1)
B = t.reshape(1, -1, 1)
T = t.reshape(1, 1, -1)
xt = x(t).reshape(1, 1, -1)
return np.sum(xt * p((T - B) / A), axis=2) / (np.abs(A.reshape(-1, 1)) ** 0.5) * dt
def icwt(xw, e, t, a):
da = a[1] - a[0]
db = t[1] - t[0]
A = a.reshape(-1, 1, 1)
B = t.reshape(1, -1, 1)
T = t.reshape(1, 1, -1)
Xw = xw.reshape(*xw.shape, 1)
return np.sum(Xw * e((T - B) / A) / (np.abs(A) ** 0.5) / (A ** 2), axis=(0, 1)) * da * db
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}$
まずはCWTを行う.これに関しては過去に作って解説もあるのでそちらを見てください.
参考: Python: ウェーブレット変換の実装：②CWT編（修正＋解説追加版）
また, cwtの関数に関しては数式まんまになるようにちょっと改良を加えていますが, やっていることは同じです.
以下CWTの実行結果
t = np.linspace(-1, 1, 200, endpoint=False)
a = np.linspace(0.01, 0.5, 100, endpoint=False)
x = lambda t: np.cos(2 * np.pi * 7 * t) + np.sin(np.pi * 6 * t)
print('x(t)')
plt.plot(t, x(t))
plt.show()
xw = cwt(x, ricker, t, a)
print('Xw(a, b)')
plt.imshow(xw, extent=[-1, 1, 0.01, 0.5], cmap='PRGn', aspect='auto',
vmax=abs(xw).max(), vmin=-abs(xw).max())
plt.show()
x(t)
Xw(a, b)
ICWT
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}$
の式で,この変換は先で行った.
次にこの式で求まったXw(a, b)から逆変換する式は
${\displaystyle x(t)=C_{\psi }^{{-1}}\int_{{-\infty }}^{{\infty}}\int_{{-\infty }}^{{\infty}}X_{w}(a,b)\frac{1}{|a|^{{1/2}}}\eta\left({\frac {t-b}{a}}\right),db\ {\frac {da}{a^{2}}}}$
でこれを用いて元の波形に逆変換する.（今回は逆変換に必要な${\eta}$にもrickerを用いる）
また,この式の
$C_{\psi } x(t)={\displaystyle \int_{{-\infty }}^{{\infty}}\int_{{-\infty }}^{{\infty}}X_{w}(a,b)\frac{1}{|a|^{{1/2}}}{\eta}\left({\frac {t-b}{a}}\right),db\ {\frac {da}{a^{2}}}}$
の部分が元の波形の形状を復元する部分で元波形はその定数倍になっています.
この定数は
${\displaystyle C_{\psi }=\int _{{-\infty }}^{{\infty }}{\frac {\overline {\hat {\psi }}(\omega ){\hat {{\eta}}}(\omega )}{|\omega |}},d\omega }$
で求められる.
※この定数はまだ求める関数ができていないので実装出来次第追記します.
また, この式が1となるような${\psi}$と${\eta}$の場合.
$x(t)={\displaystyle \int_{{-\infty }}^{{\infty}}\int_{{-\infty }}^{{\infty}}X_{w}(a,b)\frac{1}{|a|^{{1/2}}}{\eta}\left({\frac {t-b}{a}}\right),db\ {\frac {da}{a^{2}}}}$
となる.
今回実装した逆変換は上の式で変換を行うように実装している.
したがって, 実行結果を${C_{\psi }}$で割るか, ${C_{\psi }}$が1となるような${\psi}$と${\eta}$で変換する必要がある.
実験
とりあえず定数倍のまま実行してみる.
以下実行結果
plt.plot(t, icwt(xw, ricker, t, a))
plt.show()
定数を求める方法はこんなイメージ・・・だが多分正確ではない（ψとηにハットが付いてるので）
${\displaystyle C_{\psi }=\int _{{-\infty }}^{{\infty }}{\frac {\overline {\hat {\psi }}(\omega ){\hat {{\eta}}}(\omega )}{|\omega |}},d\omega }$
# np.sum(x(t) * icwt(xw, ricker, t, a) * (t[1] - t[0])) / 2
np.mean(x(t) * icwt(xw, ricker, t, a))
>>> 65.16417060673811
定数倍して元の波形と比較してみる.
plt.plot(t, icwt(xw, ricker, t, a) / 65.16417060673811)
plt.plot(t, x(t))
plt.show()
ちなみに${C_{\psi }}$が1となるような${\psi}$と${\eta}$とは一方を${C_{\psi }}$で割った関係ともいえる.
したがってこんな変換でも戻せる.
e = lambda t: ricker(t) / 65.16417060673811
plt.plot(t, icwt(xw, e, t, a))
plt.plot(t, x(t))
plt.show()
まとめ
今回は以前作成したCWTを改良し, そのICWTを実装した.
今後
定数を求める部分を改良する！
感想
ほぼほぼ理解できたので良かった(´∀｀)
追記
定数をざっくり求める処理で求めるように変更しました.
cwtの結果にdtをかけ忘れていたので修正しました.