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

　今回はPythonのPywaveletsモジュールを用いたウェーブレット変換について紹介する。

1. Pywaveletモジュール
2. 連続ウェーブレット変換
3. 離散ウェーブレット変換と逆変換
4. まとめ

1. Pywaveletモジュール

前回はScipyのcwt関数を用いた連続ウェーブレット変換について紹介したが

atatat.hatenablog.com

今回はPywaveletsモジュール

PyWavelets - Wavelet Transforms in Python — PyWavelets Documentation

を用いたウェーブレット変換を行う。Scipyとの差は、基底関数の種類、離散ウェーブレット変換、２Dウェーブレット変換などウェーブレット変換に関して多くの関数が利用できる点にある。

Scipyと比較するためウェーブレット変換を適用する信号は前回と同じとする。

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'])

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

　Pywaveletsモジュールがインストールされていれば、import pywtとすればロードできる。連続ウェーブレット変換はScipyと同名のcwt関数を使用する。引数の順番がことなるが、Scipyとほぼ同じ使い勝手でpywt.cwt('sig, width, 'mexh')とすると基底関数をメキシカンハット関数(Scipyではricker)とするウェーブレット変換を実行できる。widthについてはScipyと同じで1から30のアレイで指定する。

import pywt

width = np.arange(1, 31)
cwtmatr,freq = pywt.cwt(sig, width, 'mexh')
plt.imshow(cwtmatr,aspect='auto')
plt.colorbar()

連続ウェーブレット変換の結果もScipyと同じになる。

Pywaveletsモジュールでサポートされている基底関数を以下のように調べることができる。wavelist関数の引数kindは連続ウェーブレットではcontinuous、離散ウェーブレットではdiscreteを指定する。（両方の場合はall）

wavlist = pywt.wavelist(kind='continuous')
wavlist
>['cgau1',
 'cgau2',
 'cgau3',
 'cgau4',
 'cgau5',
 'cgau6',
 'cgau7',
 'cgau8',
 'cmor',
 'fbsp',
 'gaus1',
 'gaus2',
 'gaus3',
 'gaus4',
 'gaus5',
 'gaus6',
 'gaus7',
 'gaus8',
 'mexh',
 'morl',
 'shan']

試しにリストにあるShannon関数shanを用いてウェーブレット変換しようとすると、帯域幅と中心周波数を関数名で指定するように警告が表示される。

width = np.arange(1, 31)
cwtmatr,freq = pywt.cwt(sig, width, 'shan')
plt.imshow(np.abs(cwtmatr),aspect='auto')
plt.colorbar()
>/usr/local/lib/python3.7/dist-packages/pywt/_cwt.py:117: FutureWarning: Wavelets from the family shan, without parameters specified in the name are deprecated. The name should take the form shanB-C where B and C are floats representing the bandwidth frequency and center frequency, respectively (example: shan1.5-1.0).
  wavelet = DiscreteContinuousWavelet(wavelet)

警告の例の通りにshan1.5-1.0とすると、警告なしで実行される。

width = np.arange(1, 31)
cwtmatr,freq = pywt.cwt(sig, width, 'shan1.5-1.0')
plt.imshow(np.abs(cwtmatr),aspect='auto')
plt.colorbar()

"Shannon関数によるウェーブレット変換" — Shannon関数によるウェーブレット変換

3. 離散ウェーブレット変換と逆変換

　Pywaveletsモジュールを使用すると離散ウェーブレット変換をdwt関数を用いてdwt(sig.'db1')を実行し係数を得ることができる。ここでdb1は基底関数でどのような種類があるかはpywt.wavelist(kind='discrete')で確認できる。

width = np.arange(1, 31)
cA,cD = pywt.dwt(sig,  'db1')
plt.figure(figsize=(10,5))
plt.plot(cA)
plt.plot(cD)

フーリエ変換と同様に逆変換により元の信号に戻すことができる。離散ウェーブレット変換により得られた係数cA,cDを用いて離散ウェーブレット逆変換をidwt関数で実行すると、信号を復元することができる。

sig_r1 = pywt.idwt(cA, None ,'db1','smooth')
sig_r2 = pywt.idwt(None, cD ,'db1','smooth')
plt.figure(figsize=(10,5))
plt.plot(sig_r1+sig_r2)
plt.plot(sig,'r--')
plt.legend(['IDWT','Original'])