PyWaveletsに関するメモ (公式資料の低質の和訳)．
PyWaveletsとは
インストール
- バイナリ(pip, anaconda)
- ソースからビルド
API
Wavelet ファミリー
- Waveletファミリーは 以下のファミリー(mother wavelet(基底信号？)の大きな分類？)がbuilt-inで実装されている．
- 各ファミリーがどのような信号か知りたい人は，Matlab Wavelet Toolboxのこのページを見ると良い．
- Haar (haar)
- Daubechies (db)
- Symlets (sym)
- Coiflets (coif)
- Biorthogonal (bior)
- Reverse biorthogonal (rbio)
- “Discrete” FIR approximation of Meyer wavelet (dmey)
- Gaussian wavelets (gaus)
- Mexican hat wavelet (mexh)
- Morlet wavelet (morl)
- Complex Gaussian wavelets (cgau)
- Shannon wavelets (shan)
- Frequency B-Spline wavelets (fbsp)
- Complex Morlet wavelets (cmor)
- pywt.families()を使うと，built-inのファミリーをstrのlistとして取得できる．
>>> families = pywt.families() >>> print(families) ['haar', 'db', 'sym', 'coif', 'bior', 'rbio', 'dmey', 'gaus', 'mexh', 'morl', 'cgau', 'shan', 'fbsp', 'cmor'] >>> families_not_short = pywt.families(short=False) >>> print(families_not_short) ['Haar', 'Daubechies', 'Symlets', 'Coiflets', 'Biorthogonal', 'Reverse biorthogonal', 'Discrete Meyer (FIR Approximation)', 'Gaussian', 'Mexican hat wavelet', 'Morlet wavelet', 'Complex Gaussian wavelets', 'Shannon wavelets', 'Frequency B-Spline wavelets', 'Complex Morlet wavelets']
使用可能なWavelet
- pywt.wavelist()を使うと，指定したファミリーのwaveletをstrのlistとして取得できる．
- kind引数を'continuous', 'discrete'に指定すると，それぞれ連続，離散waveletの一覧を取得できる．
>>> wavelets_coif = pywt.wavelist('coif')
>>> print(wavelets_coif)
['coif1', 'coif2', 'coif3', 'coif4', 'coif5', 'coif6', 'coif7', 'coif8', 'coif9', 'coif10', 'coif11', 'coif12', 'coif13', 'coif14', 'coif15', 'coif16', 'coif17']
>>> wavelets_continuous = pywt.wavelist(kind='continuous')
>>> print(wavelets_continuous)
['cgau1', 'cgau2', 'cgau3', 'cgau4', 'cgau5', 'cgau6', 'cgau7', 'cgau8', 'cmor', 'fbsp', 'gaus1', 'gaus2', 'gaus3', 'gaus4', 'gaus5', 'gaus6', 'gaus7', 'gaus8', 'mexh', 'morl', 'shan']
Waveletオブジェクト
pywt.Wavelet()を使うと，離散Waveletの属性(特徴)を表すオブジェクトを取得できる．
- 連続waveletについてはpywt.ContinuousWaveletを使う．
- wavelet名の指定はpywt.wavelist()で得られる名前を使う．
- 自作waveletのWaveletオブジェクトを作るためには自身が定義したフィルタセットを filter_bank パラメータと一緒に与える必要がある．
- Waveletオブジェクトの各パラメータの意味は以下の通り
- フィルタは全てfloatのlist
- name: Wavelet名 (pywt.wavelist()で得られる名前)
- orthogonal: 直交性 (直交性を持つならTrue)
- biorthogonal: 双直交性 (双直交性を持つならTrue)
- symmetry: 対称性 (asymmetric, near symmetric, symmetricのどれか)
- 分解 (decomposition)
- dec_len: 分解フィルタ長
- dec_lo: 分解フィルタ (低周波)
- dec_hi: 分解フィルタ (高周波)
- 再構成 (reconstruction)
- rec_len: 再構成フィルタ長
- rec_lo: 再構成フィルタ (低周波)
- rec_hi: 再構成フィルタ (高周波)
- filter_bank: 全フィルタ
- listのtuple
- tupleの順番はdec_lo, dec_hi, rec_lo, rec_hi
- inverse_filter_bank: reverse wavelet filters係数のリスト(と書いてある)
- listのtuple
- tupleの順番はrec_lo[::-1], rec_hi[::-1], dec_lo[::-1], dec_hi[::-1]
- vanishing_moments_psi: wavelet関数の消失モーメントの数
- vanishing_moments_phi: scaling関数の消失モーメントの数
- 以下のDWTとCWTの取得方法が分からなかった...
>>> wavelet = pywt.Wavelet('haar')
>>> print(wavelet)
Wavelet haar
Family name: Haar
Short name: haar
Filters length: 2
Orthogonal: True
Biorthogonal: True
Symmetry: asymmetric
DWT: True
CWT: False
DWT
- PyWaveletsはWavelet変換の解像度と対象信号の次元という2つの多様性を持つ
- 解像度
- 単解像度 (dwt, idwt)
- 多重解像度 (wavedec, waverec)
- 対象信号の次元
- 1次元:
- 2次元
- n次元
- 解像度
- 上記の4個の関数名は1次元信号用
- 2次元，n次元用の関数名は末尾に2,nを付ける (例: 2次元信号の多重解像度分解は wavedec2)
1次元単解像度Wavelet変換
- dwtの返戻値は要素数が2のtuple．
- 要素0はapproximationの係数 (float64，1次元のnumpy.ndarray型)
- 要素1はdetailの係数 (float64，1次元のnumpy.ndarray型)
- idwtの引数はapproximation, detail係数とmother wavelet
>>> signal = [1, 2, 3, 4, 5, 6] >>> mother_wavelet = 'db1' >>> (cA, cD) = pywt.dwt(signal, mother_wavelet) >>> print(type(cA), type(cD)) <class 'numpy.ndarray'> <class 'numpy.ndarray'> >>> print(cA.shape, cD.shape) (3,) (3,) >>> print(cA.dtype, cD.dtype) float64 float64 >>> signal_recon = pywt.idwt(cA, cD, mother_wavelet) >>> print(signal_recon) [1. 2. 3. 4. 5. 6.] >>> print(type(signal_recon)) <class 'numpy.ndarray'> >>> print(signal_recon.shape) (6,) >>> print(signal_recon.dtype) float64
2次元単解像度Wavelet変換
- dwt2の返戻値は要素数が2のtuple．
- 要素0はapproximationの係数
- float64，2次元のnumpy.ndarray型
- 低周波成分に該当
- 要素1はdetailの係数をまとめた要素3のtuple
- tupleの要素は横方向，縦方向，斜め方向の高周波成分
- 各要素はfloat64，2次元のnumpy.ndarray型
- 要素0はapproximationの係数
- idwt2の引数はwavelet係数とmother wavelet
- wavelet係数はdwt2の返戻値と同じ形式でまとめる
>>> signal = np.random.rand(8, 8)
>>> mother_wavelet = 'db1'
>>> c = pywt.dwt2(signal, mother_wavelet)
>>> print('c:', type(c), len(c))
c: <class 'tuple'> 2
>>> print(' c[0]:', type(c[0]), c[0].shape)
c[0]: <class 'numpy.ndarray'> (4, 4)
>>> print(' c[1]:', type(c[1]), len(c[1]))
c[1]: <class 'tuple'> 3
>>> [print(' c[1][%d]:' % n, type(c[1][n]), c[1][n].shape, c[1][n].dtype) for n in range(len(c[1]))]
c[1][0]: <class 'numpy.ndarray'> (4, 4) float64
c[1][1]: <class 'numpy.ndarray'> (4, 4) float64
c[1][2]: <class 'numpy.ndarray'> (4, 4) float64
>>> cA, (cH, cV, cD) = c
>>> signal_recon = pywt.idwt2((cA, (cH, cV, cD)), mother_wavelet)
>>> print(np.testing.assert_almost_equal(signal, signal_recon))
None
>>> print(type(signal_recon))
<class 'numpy.ndarray'>
>>> print(signal_recon.shape)
(8, 8)
>>> print(signal_recon.dtype)
float64
多重解像度分解
- wavedec2の返戻値は要素数が引数levelで指定した解像度と同数のlist．
- 要素0はapproximationの係数
- float64，2次元のnumpy.ndarray型
- 低周波成分に該当
- 要素1〜level-1はdetailの係数をまとめた要素3のtuple
- tupleの要素は横方向，縦方向，斜め方向の高周波成分
- 各要素はfloat64，2次元のnumpy.ndarray型
- 要素0はapproximationの係数
- waverec2は調べてないけど，idwt2と似た感じだと思う
level=2とした時の画像とwavelet係数の対応関係を図にしてみた(分かりにくい)．
[c[0], (c[1][0], c[1][1], c[1][2]), (c[2][0], c[2][1], c[2][2])] <---> ------------------------------------------ | | | | | c[0] | c[1][0] | | | | | | --------------------- c[2][0] | | | | | | c[1][1] | c[1][2] | | | | | | ------------------------------------------ | | | | | | | | | | c[2][1] | c[2][2] | | | | | | | | | | ------------------------------------------
>>> image = pywt.data.camera().astype(np.float32)
>>> print(' ', type(image), image.shape)
<class 'numpy.ndarray'> (512, 512)
>>> mother_wavelet = 'db1'
>>> for level in range(5):
>>> c = pywt.wavedec2(image, mother_wavelet, mode='periodization', level=level)
>>> print('Wavelet info with level=%d' % level)
>>> print(' level=%d' % level, type(c), len(c))
level=0 <class 'list'> 1
level=1 <class 'list'> 2
level=2 <class 'list'> 3
level=3 <class 'list'> 4
level=4 <class 'list'> 5
>>> level = 2
>>> c = pywt.wavedec2(image, mother_wavelet, mode='periodization', level=level)
>>> print('level=%d' % level)
>>> print(' c[0]:', type(c[0]), c[0].shape)
level=2
c[0]: <class 'numpy.ndarray'> (128, 128)
>>> for l in range(level):
>>> print(' c[%d]' % (l+1), type(c[l+1]), len(c[l+1]))
>>> for elem in range(len(c[l+1])):
>>> print(' c[%d][%d]' % (l+1, elem), type(c[l+1][elem]), c[l+1][elem].shape)
c[1] <class 'tuple'> 3
c[1][0] <class 'numpy.ndarray'> (128, 128)
c[1][1] <class 'numpy.ndarray'> (128, 128)
c[1][2] <class 'numpy.ndarray'> (128, 128)
c[2] <class 'tuple'> 3
c[2][0] <class 'numpy.ndarray'> (256, 256)
c[2][1] <class 'numpy.ndarray'> (256, 256)
c[2][2] <class 'numpy.ndarray'> (256, 256)
Wavelet係数の扱い方
- Wavelet係数と単一の(元信号と同じ次元の)numpy.ndarrayの相互変換
- 多重解像度のWavelet係数を扱いやすくする方法
- pywt.coeffs_to_array: Wavelet係数からarrayへ変換
- coeff_arr (第1返戻値): 原信号と同じ次元数のnumpy.ndarray型オブジェクト．
- coeff_slices (第2返戻値): 要素数がlevel引数で指定した多重解像度と同じlist型オブジェクト．
- pywt.array_to_coeffs: arrayからWavelet係数へ変換
- pywt.coeffs_to_array: Wavelet係数からarrayへ変換
- pywt.coeffs_to_arrayによってwavelet係数を1枚の画像に変換できる
- levelを変えてもpywt.coeffs_to_arrayによって得られる画像の解像度は変わらない (画像サイズが2の階乗でない場合は違うかも)
- levelの数が多重解像度の数であることが分かる
>>> image = pywt.data.camera().astype(np.float32)
>>> print('Image info: , type(image), image.shape)
Image info: <class 'numpy.ndarray'> (512, 512)
>>> mother_wavelet = 'db1'
>>> levels = 4
>>> fig, axes = plt.subplots(1, levels+1, figsize=[12, 10])
>>> axes[0].imshow(image, cmap=plt.cm.gray)
>>> axes[0].set_axis_off()
>>> axes[0].set_title('Original image')
>>> for level in range(4):
>>> print('level=%d' % level)
>>> coeff = pywt.wavedec2(image, mother_wavelet, mode='periodization', level=level)
>>> arr, coeff_slices = pywt.coeffs_to_array(coeff)
>>> print(' coeff:', type(arr), arr.shape)
>>> print(' array:', type(coeff_slices), len(coeff_slices))
>>> axes[level+1].imshow(arr, cmap=plt.cm.gray)
>>> axes[level+1].set_axis_off()
>>> axes[level+1].set_title('Wavelet coeff. at level %d' % level)
level=0
coeff: <class 'numpy.ndarray'> (512, 512)
array: <class 'list'> 1
level=1
coeff: <class 'numpy.ndarray'> (512, 512)
array: <class 'list'> 2
level=2
coeff: <class 'numpy.ndarray'> (512, 512)
array: <class 'list'> 3
level=3
coeff: <class 'numpy.ndarray'> (512, 512)
array: <class 'list'> 4
- Wavelet係数をまとめた画像(信号)から各解像度，方向の成分を取り出すには，pywt.coeffs_to_arrayの第2返戻値を使う
- coeff_slices[level][key]
- 解像度level，信号の方向keyのWavelet係数をcoeffから取得するためのslice型オブジェクトを要素に持つtuple型オブジェクト．
- keyは'da', 'ad', 'dd'の3種類でそれぞれ横，縦，斜め方向の係数を表す
- coeff_slices[level][key]
>>> image = pywt.data.camera().astype(np.float32)
>>> mother_wavelet = 'db1'
>>> level = 2
>>> coeff = pywt.wavedec2(image, mother_wavelet, mode='periodization', level=level)
>>> arr, coeff_slices = pywt.coeffs_to_array(coeff)
>>> fig, axes = plt.subplots(level+2, 3, figsize=[12, 10])
>>> axes[0, 0].imshow(image, cmap=plt.cm.gray)
>>> axes[0, 0].set_axis_off()
>>> axes[0, 0].set_title('Original image')
>>> axes[0, 1].imshow(arr, cmap=plt.cm.gray)
>>> axes[0, 1].set_axis_off()
>>> axes[0, 1].set_title('Wavelet coefficients')
>>> axes[0, 2].set_axis_off()
>>> print(' coeff:', type(arr), arr.shape)
>>> print(' array:', type(coeff_slices), len(coeff_slices))
coeff: <class 'numpy.ndarray'> (512, 512)
array: <class 'list'> 3
>>> for l in range(level):
>>> print(' array for level=%d info' % (l + 1), type(coeff_slices[l + 1]), len(coeff_slices[l + 1]))
>>> if l == 0:
>>> axes[1, 0].imshow(arr[coeff_slices[l]], cmap=plt.cm.gray)
>>> axes[1, 0].set_title('Low frequency')
>>> axes[1, 0].set_axis_off()
>>> axes[1, 1].set_axis_off()
>>> axes[1, 2].set_axis_off()
>>> for n, key in enumerate(coeff_slices[l+1].keys()):
>>> axes[l+2, n].imshow(arr[coeff_slices[l+1][key]], cmap=plt.cm.gray)
>>> axes[l+2, n].set_title('%s at level %d' % (key, l))
>>> axes[l+2, n].set_axis_off()
>>> print(' ', key, type(coeff_slices[l+1][key]), coeff_slices[l+1][key])
>>> plt.show()
array for level=1 info <class 'dict'> 3
ad <class 'tuple'> (slice(None, 128, None), slice(128, 256, None))
da <class 'tuple'> (slice(128, 256, None), slice(None, 128, None))
dd <class 'tuple'> (slice(128, 256, None), slice(128, 256, None))
array for level=2 info <class 'dict'> 3
ad <class 'tuple'> (slice(None, 256, None), slice(256, 512, None))
da <class 'tuple'> (slice(256, 512, None), slice(None, 256, None))
dd <class 'tuple'> (slice(256, 512, None), slice(256, 512, None))
