今回はPythonのPywaveletsモジュールを用いたウェーブレット変換について紹介する。
1. Pywaveletモジュール
前回はScipyのcwt関数を用いた連続ウェーブレット変換について紹介したが
今回はPywaveletsモジュール
PyWavelets - Wavelet Transforms in Python — PyWavelets Documentation
を用いたウェーブレット変換を行う。Scipyとの差は、基底関数の種類、離散ウェーブレット変換、2Dウェーブレット変換などウェーブレット変換に関して多くの関数が利用できる点にある。
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()
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'])
4. まとめ
今回はPywaveletsモジュールを用いた連続ウェーブレット変換、離散ウェーブレット変換について紹介した。