Pythonによるデータ処理9 ～ウェーブレット変換

　今回はPythonによるウェーブレット変換について紹介する。

1. ウェーブレット変換
2. 信号生成
3. 連続ウェーブレット変換
3. まとめ

1. ウェーブレット変換

　ウェーブレット変換はフーリエ変換のような周波数解析の１種で基底関数としてウェーブレット関数Ψ (t)を用いて、連続ウェーブレット変換は以下のように積分変換する。

$\displaystyle X (a,b) = \frac{1}{a} \int^{\infty}_{-\infty} x(t) \Psi^{*} (\frac{t-b}{a} ) dt$

基本的には信号がウェーブレット関数の拡大・縮小、平行移動によってあらわす。ウェーブレット関数は許容条件がありメキシカンハット関数などが使われる。

2. 信号生成

ウェーブレット変換の例として、以下のサンプルコードを基にテスト用の信号を作成する。

scipy.signal.cwt — SciPy v1.7.1 Manual

テスト用信号sigは、cos波tri、ガウシアンパルスpulse、白色ノイズnoiseを加算したものとする。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
%matplotlib inline

t = np.linspace(-1, 1, 200, endpoint=False)
tri = np.cos(2 * np.pi * 7 * t)
pulse = signal.gausspulse(t - 0.4, fc=2)
noise = np.random.normal(0,0.5,np.size(t))

sig  =  tri+ pulse +noise 

plt.figure(figsize=(8,6))
plt.plot(t,sig)
plt.plot(t,tri-3,'r')
plt.plot(t,pulse-6,'b')
plt.plot(t,noise-8,'g')
plt.legend(['sig','tri','pulse','noise'])

上記のコードを実行して、テスト用の信号、cos波、ガウシアンパルス、ノイズをそれぞれ表示した。

3. 連続ウェーブレット変換

　Pythonのウェーブレット変換できるモジュールはいくつかあるが、今回はScipyのcwt関数を利用する。使い方としては上記で生成したテスト用信号sigとしてcwtmatr = signal.cwt(sig, signal.ricker, width)とするとウェーブレット関数としてricker (メキシカンハット関数)、スケーリングパラメータ:widthで連続ウェーブレット変換を実行する。widthはアレイで生成し、np.arange(1,31)として1から30までの値としている。

widths = np.arange(1, 31)
cwtmatr = signal.cwt(sig, signal.ricker, widths)
plt.imshow(cwtmatr,aspect='auto')
plt.colorbar()

cwt関数を実行して、imshowを用いてマップ表示する。スケーリング係数が5ぐらいまではテスト用信号に含まれるcos波成分、スケーリング係数10ぐらいはガウシアンパルス成分をよくあらわしている。

Scipyのウェーブレット基底用関数としてricker以外にmorletも使用できる。ただし複素数で出力されるので、実数に変換して表示する必要がある。

cwtmatr = signal.cwt(sig, signal.morlet, widths)
plt.imshow(np.abs(cwtmatr),aspect='auto')
plt.colorbar()

"ウェーブレット変換, morlet関数" — ウェーブレット変換, morlet関数

ウェーブレット変換に使用したRicker関数を直接プロットするとスケーリング係数と関数の形の関係がわかる。ricker(points, a)のpointsはプロット点数、aはスケーリングパラメータになり、a=1だと幅が狭、a=30だと幅が広いメキシカンハット関数になっていることが確認できる。

plt.plot(signal.ricker(100,1),'g')
plt.plot(signal.ricker(100,10),'r')
plt.plot(signal.ricker(100,30),'b')

Ricker関数と同じようにプロット点数とスケーリングパラメータを引数とするユーザ定義関数をウェーブレット変換に使用することができる。以下は指数関数をウェーブレット変換に使用したが、指数関数自体は基底関数の許容条件を満たしておらず、単にユーザ定義関数がcwt関数によるウェーブレット変換に使えるということのみを示している。　　　

def mexp(points,a):
  x = np.linspace(-1,1,points)
  y = np.exp(-x**2/(a/50)**2)
  return y

plt.figure(figsize=(12,4))

plt.subplot(1,2,1)
cwtmatr = signal.cwt(sig, mexp, widths)
plt.imshow(cwtmatr,aspect='auto')
plt.colorbar()
plt.subplot(1,2,2)
plt.plot(mexp(200,1),'g')
plt.plot(mexp(200,10),'r')
plt.plot(mexp(200,30),'b')