Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

はじめに

目標

本投稿では, 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変換後である.

サンプル版.png

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)

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

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

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

大まかなイメージ.png

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

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

自作版.png

まとめ

scpyのcwtと(ほぼ)同じCWTを実装した.

感想

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

sai-sui
いつもはDeep Learning関連のこと色々やってる よく使う(必要になる)波形, 画像処理系を勉強していきたい
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away