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

はじめに

今回はウェーブレット変換の実装：②CWT編で作ったCWTの改良を行う.
改良点
- ②で作成したrikcer関数をそのまま引数で与えてMother Waveletに設定可能にする.
- Mother Waveletの基本サイズを引数rで調整できるようにする.
- sigmoid関数を追加する.（同様に変換できるかも実験）

改良したコード

sig.py

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 sigmoid(t, s=1.0):
    return (np.tanh(s * t / 2) + 1) / 2


def cwt(y, mw, S, r=10.0):
    out = []
    n = len(y) / (2 * r)
    t = np.linspace(-n, n, len(y), endpoint=False)
    for s in S:
        wave = mw(t, s)
        out.append(np.convolve(y, wave, mode='same'))
    return np.array(out)

実験（遊んでみる）

準備

上のsig.pyをとってきてはじめのおまじないを唱える(∩｡･ｏ･｡)っ.ﾟ☆｡'`

import matplotlib.pyplot as plt
import numpy as np

from sig import *

適当な合成波を作る

ウェーブレット変換の実装：②CWT編で作ったのと同じやつ.

t = np.linspace(-1, 1, 200, endpoint=False)
y  = np.cos(2 * np.pi * 7 * t) + np.sin(np.pi * 6 * t)
plt.plot(t, y)
plt.show()

ricker関数の引数sを変更しながら色んな波を作り表示してみる.

t = np.linspace(-5, 5, 100, endpoint=False)
S = np.linspace(0.1, 1, 10, endpoint=False)
for s in S:
    plt.plot(t, ricker(t, s))
plt.show()

上のrickerをMother Waveletに設定してCWTを実行.

sp = 0.1
ep = 3.0
reso = 100
r = 10
mw = ricker
widths = np.linspace(sp, ep, reso, endpoint=False)
cy = cwt(y, mw, widths, r)
print('output')
plt.imshow(cy, extent=[-1, 1, ep*r, sp*r], cmap='PRGn', aspect='auto',
           vmax=abs(cy).max(), vmin=-abs(cy).max())
plt.show()

違いはあるがウェーブレット変換の実装：②CWT編で得られた結果とほぼ同じ.

続いてsigmoid関数の引数sを変更しながら色んな波を作り表示してみる.

t = np.linspace(-10, 10, 100, endpoint=False)
S = np.linspace(-1, 1, 10, endpoint=False)
for s in S:
    plt.plot(t, sigmoid(t, s))
plt.show()

rickerと同じように変換する

sp = 0.1
ep = 3.0
reso = 100
r = 10
widths = np.linspace(sp, ep, reso, endpoint=False)
cy = cwt(y, sigmoid, widths, r)
print('output')
plt.imshow(cy, extent=[-1, 1, ep*r, sp*r], cmap='PRGn', aspect='auto',
           vmax=abs(cy).max(), vmin=-abs(cy).max())
plt.show()

出来た！ρ(・ω・、)

まとめ

とりあえず目標の改良は出来た. あとは色んな関数を突っ込めば色んな変換ができる.