ウェーブレット変換とは

wavelet transform

時系列データに対して、信号の時間と周波数の関係を同時に解析するための手法。

同じく時系列データを解析する手法としてフーリエ変換（cf. fft）があるが、フーリエ変換によるスペクトル解析では、時間方向の情報が失われてしまうため、例えば「時間が後ろになるほど周波数が高くなる」ような状況を解析することが難しい。

ウェーブレット変換ではその欠点を補い、時間とともにデータの周期性が変化していく様子を捉えることができる。

基本的な考え方・変換手順

ウェーブレットの生成

信号の周波数成分を解析するための物差しとして導入する「波の断片」を ウェーブレットと呼ぶ。
基準となるマザーウェーブレット（ウェーブレット基底関数）をメキシカンハット関数

$$\psi(t) = \left(1 - t^2 \right) \exp{ \left( - \cfrac{t^2}{2} \right)}$$

などで定義する（以後の説明では、ウェーブレットとしてメキシカンハット関数を使用）。

これを時間軸方向に拡大・縮小 / 平行移動して色々な形のウェーブレットを生成し、後の変換に用いる：

$$\psi_{\sigma, t_0}(t) = \cfrac{1}{\sqrt{\sigma}} \psi \left( \cfrac{t-t_0}{\sigma} \right) = \cfrac{1}{\sqrt{\sigma}} \left(1 - \cfrac{(t-t_0)^2}{\sigma^2} \right) \exp{ \left( - \cfrac{(t-t_0)^2}{2\sigma^2} \right)}$$

パラメータ	説明
$\sigma$	拡大縮小のパラメータ。波の波長に相当(逆数 $\frac{1}{\sigma}$ は周波数に相当)
$t_0$	平行移動のパラメータ。波の中心位置に相当

ここで、$\sigma$ 倍に引き伸ばされた波のエネルギーのスケールを合わせるため、重み $\frac{1}{\sqrt{\sigma}}$ をかけている。

（展開：作図に使ったコード）

	from matplotlib import pyplot as plt
	import numpy as np
	def mexican_hat(t, sigma, t0):
	tmp = ((t - t0) / sigma)**2
	return (1 - tmp) * np.exp(-tmp/2) / np.sqrt(sigma)
	t = np.arange(2000) / 100
	t0_list = [5.0, 10.0, 15.0]
	sigma_list = [0.5, 1.0, 2.0]
	fig = plt.figure(figsize=(9,9))
	c = 0
	for t0 in t0_list:
	for sigma in sigma_list:
	c += 1
	y = mexican_hat(t, sigma, t0)
	ax = fig.add_subplot(3, 3, c)
	ax.set_title('$\sigma = {}, t_0 = {}$'.format(sigma, t0))
	ax.set_ylim(-1.0, 1.5)
	ax.plot(t, y)
	ax.grid()
	plt.tight_layout()
	plt.show()

view raw wavelet.py hosted with ❤ by GitHub

時間・周波数的な特徴の計算

用意した色々なウェーブレット $\psi_{\sigma, t_0}(t)$ と元の時系列データ $f(t)$ とで、同じ時刻の成分どうしの積 $\psi_{\sigma, t_0}(t) f(t)$ を取り、全時間で積分することを考える。

$$I(\sigma, t_0) = \int_{-\infty}^{\infty} \psi_{\sigma, t_0}(t) f(t) dt$$

ウェーブレット $\psi_{\sigma, t_0}(t)$ は、以下の特徴を持つ波であると考えることができる。

• 特定の時刻 $t_0$ 付近だけに高い振幅を持つ
• 波長（周期）が $\sigma$ に比例する

なので、

• 積 $\psi_{\sigma, t_0}(t) f(t)$ の値は、時刻 $t_0$ 付近以外ではゼロに近い値になる
• その積分値 $I(\sigma, t_0)$ は、時刻 $t_0$ 付近に同じ波長・同じ位相の波が存在すると大きな値になる
  • 位相が同じでも、ウェーブレットと時系列データの波長が異なると、積分値は打ち消しあって小さくなる

時系列データに対して様々なウェーブレットをかけて積分した具体例を下図に示す。

• blue：解析対象の時系列データ $y = \sin(t^2/5) \ \ \ (0 \le t \le 20)$
• orange：ウェーブレット $y = \psi_{\sigma, t_0}(t)$
• green：時系列データとウェーブレットの積関数

（展開：作図に使ったコード）

	from matplotlib import pyplot as plt
	import numpy as np
	def mexican_hat(t, sigma, t0):
	tmp = ((t - t0) / sigma)**2
	return (1 - tmp) * np.exp(-tmp/2) / np.sqrt(sigma)
	t = np.arange(2000) / 100
	y = np.sin(t**2 / 5)
	t0_list = [2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 17.5]
	sigma_list = [0.25, 0.5, 0.75, 1.0, 1.25, 1.5]
	fig = plt.figure(figsize=(12,9))
	c = 0
	for sigma in sigma_list:
	for t0 in t0_list:
	c += 1
	wavelet = mexican_hat(t, sigma, t0)
	ax = fig.add_subplot(len(sigma_list), len(t0_list), c)
	ax.set_ylim(-1.5, 1.5)
	ax.axes.xaxis.set_visible(False)
	ax.axes.yaxis.set_visible(False)
	ax.plot(t, y, lw=0.5)
	ax.plot(t, wavelet, lw=0.7)
	ax.plot(t, wavelet * y)
	ax.grid()
	ax.text(0, 1.2, '$\sigma = {}, t_0 = {}$'.format(sigma, t0), fontsize=8)
	I = (wavelet*y).sum()
	ax.text(0, 0.8, '$I = {:.3f}$'.format(I), fontsize=10)
	plt.tight_layout()
	plt.show()

view raw wavelet-x-timeseries.py hosted with ❤ by GitHub

ウェーブレット変換

もっと細かく $\sigma, t_0$ を区切って積分 $I(\sigma, t_0)$ を計算し、2次元ヒートマップを作成することで、時間とともに波長（周波数）が変化していく様子を捉えることができる。

	from matplotlib import pyplot as plt
	import numpy as np
	def mexican_hat(t, sigma, t0):
	tmp = ((t - t0) / sigma)**2
	return (1 - tmp) * np.exp(-tmp/2) / np.sqrt(sigma)
	t = np.arange(1000) / 50
	ys = [
	np.sin(t**2/7),
	np.sin(t**2/7) + np.sin(t**3/20),
	np.sin(t**2/7) * np.sin(t**3/20)
	]
	y_titles = [
	'$y = \sin{(t^2/7)}$',
	'$y = \sin{(t^2/7)} + \sin{(t^3/20)}$',
	'$y = \sin{(t^2/7)} \sin{(t^3/20)}$'
	]
	frequencies = np.linspace(0.025, 10, 400)
	fig = plt.figure(figsize=(12,8))
	for y, y_title, i in zip(ys, y_titles, range(1, 3+1)):
	Is = []
	for freq in frequencies:
	sigma = 1.0 / freq
	for t0 in t:
	wavelet = mexican_hat(t, sigma, t0)
	I = (y * wavelet).sum()
	Is.append(I)
	Is = np.array(Is).reshape(len(frequencies), len(t))
	ax1 = fig.add_subplot(2, 3, i)
	ax1.set_title(y_title)
	ax1.plot(t, y)
	ax1.set_xlabel('$t$')
	ax1.set_ylabel('$y$')
	ax2 = fig.add_subplot(2, 3, i+3)
	ax2.set_title('$I(\sigma, t_0)$')
	ax2.set_xlabel('Time $t_0$')
	ax2.set_ylabel('Frequency $1 / \sigma$')
	im = ax2.imshow(Is, cmap="autumn", extent=(t[0], t[-1], frequencies[-1], frequencies[0]), aspect=2.0)
	fig.colorbar(im, orientation='horizontal', label='color')
	plt.tight_layout()
	plt.show()

view raw wavelet-transform.py hosted with ❤ by GitHub