書記が数学やるだけ#820 ウェーブレット解析

ウェーブレット解析に必要な計算などをまとめておく。

問題

説明

フーリエ変換では専ら周波数特性が解析され，時間領域の情報は残らなかった。これに窓関数をかける短時間フーリエ変換では，時間領域の解析が可能になるが，窓幅を周波数に合わせて固定する必要があるため，広い周波数領域の解析には向かなかった。これを解決するのがウェーブレット変換で，基底関数の拡大縮小を行うので，広い周波数領域の解析が可能となる。

ウェーブレットとは波を拡大縮小・平行移動して足し合わせたものであり，これを用いて変換する。

離散ウェーブレット変換について，特にハールウェーブレットが用いられる。

解答

ウェーブレット級数展開を以下に示しておく。

ハールウェーブレットの概形を以下に示す，元の形に対し拡大縮小・平行移動を行うことができる。

ハールウェーブレットの正規直交性について，内積の計算から示す。

まずp>=2，つまりm≠kの場合でS=0であることを示す。

次にqの値によりSの値を決める。

ハールウェーブレットによるウェーブレット級数展開の係数を計算する。

スカログラムとは信号の連続ウェーブレット変換 (CWT) の絶対値であり，以下のコードで可視化できる。

import numpy as np
import matplotlib.pyplot as plt
import pywt
import scipy

# 与えられた信号
t_start = 0
t_end = 2
num_points = 1000
tm = np.linspace(t_start, t_end, num_points)
amp = np.abs(np.sin(2.0 * np.pi * 2.0 * tm)) + np.tanh(tm)
freq = 10.0 + 10.0 * tm * tm + np.cos(2.0 * np.pi * 3.0 * tm)
sig = amp * np.sin(2.0 * np.pi * freq * tm) + amp * 5.0e-2 * np.random.randn(len(tm))

# スケール毎の連続ウェーブレット変換
wavelet = 'cmor'  # 連続モロー ウェーブレットを使用
scales = np.arange(1, 128)  # 使用するスケールの範囲
coefficients, frequencies = pywt.cwt(sig, scales, wavelet, sampling_period=1/(t_end - t_start))

# スカログラムの表示
plt.imshow(np.abs(coefficients), aspect='auto', extent=[t_start, t_end, scales[-1], scales[0]], cmap='jet', interpolation='bilinear')
plt.colorbar(label='Magnitude')
plt.title('Wavelet Transform Spectrogram')
plt.xlabel('Time')
plt.ylabel('Scale')
plt.show()

本記事のもくじはこちら：