【時間-周波数解析の基礎】Pythonでウェーブレット変換を実装【Wavelet】

【時間-周波数解析の基礎】Pythonでウェーブレット変換を実装【Wavelet】

短時間フーリエ変換 (STFT) を用いて時間周波数解析を行う場合,解析したい信号の時間変化に対して窓関数幅を適切に選択する必要がありました。
しかしながら,窓関数幅の選択は試行錯誤が伴い,手間がかかる場合があります。

これを解決する一つの方策がウェーブレット解析です。

この記事ではウェーブレット変換を用いた時間周波数解析,スカログラムの表示方法について紹介します。
なお,本記事は C. Torrence and G. P. Compo (1998) の論文に大きく依拠しています。

動作検証済み環境
macOS Big Sur, python3.7.9

連続ウェーブレット変換

連続ウェーブレット変換

ウェーブレット $\psi$ に関する,信号 $x(t)$ の連続ウェーブレット変換 (continuous wavelet transform: CWT) は,

$$
\mathcal{W}_\psi[x](s, \tau) := \int_{-\infty}^{\infty} x(t)\psi^*\left(\frac{\tau – t}{s}\right)\mathrm{d}t
$$

と定義されます。

$s$ はスケールというパラメータで周期 (したがって,その逆数の周波数) と関連します。$\tau$ はシフトというパラメータで,時刻に対応します。

フーリエ変換のときと同じように,$|\mathcal{W}_\psi[x](s, \tau)|$ は時刻 $t = \tau$のときのスケール $s$ を持つ成分の振幅であると解釈することができます。

時刻 $t = 0$ 付近に局在するスケールの異なるウェーブレット (接尾辞レット (let) は「小さいもの」を表すもの (例えば book (本) に対する booklet (小冊子)) で,波の小破片みたいなニュアンスがあると考えられます) をずらしながら掛けて積分することによって,ウェーブレットに近いものがもとの信号にある場合にはウェーブレットの値として0でないものが出てくるということになります。

ウェーブレット $\psi$ は,窓関数×正弦波というような形をしているものを用います。スケールを変更することは,周波数を変えると同時に窓関数幅を変えていることになります。

すなわち,振動数の高いものに対しては,窓関数の幅を小さくし,振動数の小さいものに対しては,窓関数の幅を大きくするという操作が自動的に行われていることになります。

STFTのスペクトログラムに対応して,横軸に時間,縦軸に周波数 (スケール) をとって表示したものをスカログラム (scalogram) といいます。

連続ウェーブレット変換の離散化

上で定義した連続ウェーブレット変換を離散的に行います。

ウェーブレットのもとになるマザーウェーブレット (mother wavelet) と呼ばれる時間の関数 $\psi_0$ を用意しておいて,

$$
W_n(s) = \sum_{n’=0}^{N-1} x_n\psi_0^*\left(\frac{(n’-n)\Delta t}{s}\right)
$$

連続ウェーブレット変換の近似計算 (離散化) を計算するためには,すべてのスケールについて $N$ 回の畳込み計算を行う必要があります。
この場合,入力信号長が長くなってしまうと計算時間が非常長くなってしまいます。

これを回避するために,畳み込み定理を使います。畳み込み定理とは,時間空間での畳み込みの計算は周波数空間での掛け算になるというものです。
畳み込み定理を利用した方法は周波数空間への変換が必要になり一手間かかるように思えますが,これにFFTアルゴリズムを使うと非常に効率よく計算でき,トータルでみると計算量を大幅に節約できます。

マザーウェーブレット

よく使用されるマザーウェーブレットとしては,Morletのウェーブレットがあります。
Morletはウェーブレット解析を考案したフランスの技術士で日本語では「モルレー」や「モレー」と表記されることが多いです。

$$
\psi_\mathrm{M}(t) = \pi^{-1/4}e^{i\omega_0t}e^{-t^2/2}
$$

と定義され,

という形の関数です。赤の実線は実部,青の破線は虚部を表しています。
後で必要になるスケール $s$ のMorletのウェーブレット$\psi_\mathrm{M}(t/s)$ のフーリエ変換は

$$
\Psi_\mathrm{M}(s\omega) = \pi^{-1/4} H(\omega)e^{-(s\omega – \omega_0)^2/2}
$$

です。

マザーウェーブレットとして用いられる関数は他にもあります。詳しくは C. Torrence and G. P. Compo (1998) の論文を参照してください。

スケールと周波数の対応

ウェーブレット解析においては,周波数に対応するスケール (scale) という概念が出てきます。これは,マザーウェーブレットをどれくらい時間軸方向に伸縮させるかというパラメータです。

マザーウェーブレットとしてMorletのウェーブレットを使用するとき,スケール $s$ と (疑似) 周波数 $f$ の対応関係は,

$$
f = \frac{\omega_0 + \sqrt{2 + \omega_0^2}}{4\pi s}
$$

となります。$\omega_0 = 6$ のとき,$f \simeq 0.968 s^{-1}$ となります。

スケールの離散化

スケールを離散化する必要があります。次のように離散化されることが多いようです:

$$
s_j = s_0\cdot 2^{j\Delta j}, \quad j = 0, 1, \dots, J , \\ J = \frac{1}{\Delta j}\log_2\left(\frac{N\Delta t}{s_0}\right)
$$

ここで,$s_0$ は解像可能な最小のスケール (ナイキスト周波数に対応するもの) で,フーリエ変換したときとの対応から,$s_0 = 2\Delta t$ とします。
$J$ は最も大きなスケールに対応するものです。
$\Delta j$ は周波数の解像度 (分割数) を決めるパラメータです。
スケール方向に十分な解像度が得られるように小さくする必要があります。

Cone of Interference (COI)

Cone of Interference は計算されたスカログラムの信頼性の低い領域を示すものです。
ウェーブレットがエッジの影響を受けてしまう場合は,その時間帯のスカログラムは信憑性が低いと考えられます。

モルレーのウェーブレットを使った場合,スケール $s$ のとき,エッジから $\sqrt{2}s$ だけ離れるとエッジの影響がないとされます。($\sqrt{2}s$ はウェーブレットの包絡線が $e^{-2}$ 程度減少する時間)。

スカログラムを作成する方法

上で説明したものをpythonでコーディングすると以下のようになります。

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patheffects as path_effects
import numpy as np

#------------------------------------------------------------------------------
# テスト信号の生成
#------------------------------------------------------------------------------
dt = 1.0e-4
tms = 0.0
tme = 2.0
tm = np.arange(tms, tme, dt)

amp = 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))

#------------------------------------------------------------------------------
# 諸々の設定
freq_upper = 300.0 # in kHz
out_fig_path = "scalogram.png" 

#------------------------------------------------------------------------------
# 連続ウェーブレット変換
# --- n_cwt をlen(sig)よりも大きな2のべき乗の数になるように設定
n_cwt = int(2**(np.ceil(np.log2(len(sig)))))

# --- 後で使うパラメータを定義
dj = 0.125
omega0 = 6.0
s0 = 2.0*dt
J = int(np.log2(n_cwt*dt/s0)/dj)

# --- スケール
s = s0*2.0**(dj*np.arange(0, J+1, 1))

# --- n_cwt個のデータになるようにゼロパディングをして,DC成分を除く
x = np.zeros(n_cwt)
x[0:len(sig)] = sig[0:len(sig)] - np.mean(sig)

# --- omega array
omega = 2.0*np.pi*np.fft.fftfreq(n_cwt, dt)

# --- FFTを使って離散ウェーブレット変換する
X = np.fft.fft(x)
cwt = np.zeros((J+1, n_cwt), dtype=complex) # CWT array 

Hev = np.array(omega > 0.0)
for j in range(J+1):
    Psi = np.sqrt(2.0*np.pi*s[j]/dt)*np.pi**(-0.25)*np.exp(-(s[j]*omega-omega0)**2/2.0)*Hev
    cwt[j, :] = np.fft.ifft(X*np.conjugate(Psi))

s_to_f = (omega0 + np.sqrt(2 + omega0**2)) / (4.0*np.pi) 
freq_cwt = s_to_f / s
cwt = cwt[:, 0:len(sig)]

# --- cone of interference
COI = np.zeros_like(tm)
COI[0] = 0.5/dt
COI[1:len(tm)//2] = np.sqrt(2)*s_to_f/tm[1:len(tm)//2]
COI[len(tm)//2:-1] = np.sqrt(2)*s_to_f/(tm[-1]-tm[len(tm)//2:-1])
COI[-1] = 0.5/dt

#------------------------------------------------------------------------------
# 解析結果の可視化
figsize = (210/25.4, 294/25.4)
dpi = 200
fig = plt.figure(figsize=figsize, dpi=dpi)

# --- 図の設定 (全体)
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['xtick.top'] = True
plt.rcParams['xtick.major.size'] = 6
plt.rcParams['xtick.minor.size'] = 3
plt.rcParams['xtick.minor.visible'] = True
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['ytick.right'] = True
plt.rcParams['ytick.major.size'] = 6
plt.rcParams['ytick.minor.size'] = 3
plt.rcParams['ytick.minor.visible'] = True
plt.rcParams["font.size"] = 14
plt.rcParams['font.family'] = 'Helvetica'

# プロット枠 (axes) の設定
ax1 = fig.add_axes([0.15, 0.55, 0.70, 0.3])
ax_sc = fig.add_axes([0.15, 0.2, 0.70, 0.30])
cx_sc = fig.add_axes([0.87, 0.2, 0.02, 0.30])

# 元データのプロット
ax1.set_xlim(tms, tme)
ax1.set_xlabel('')
ax1.tick_params(labelbottom=False)
ax1.set_ylabel('x')

ax1.plot(tm, sig, c='black')

# スペクトログラムのプロット
ax_sc.set_xlim(tms, tme)
ax_sc.set_xlabel('time (s)')
ax_sc.tick_params(labelbottom=True)
ax_sc.set_ylim(0, freq_upper)
ax_sc.set_ylabel('frequency\n(Hz)')

# カラーバーのレンジの設定
norm = mpl.colors.Normalize(vmin=np.log10(np.abs(cwt[freq_cwt < freq_upper, :])**2).min(), 
                            vmax=np.log10(np.abs(cwt[freq_cwt < freq_upper, :])**2).max())
# カラーマップ
cmap = mpl.cm.jet

# スカログラムのプロット
ax_sc.contourf(tm, freq_cwt, np.log10(np.abs(cwt)**2), 
                norm=norm, levels=256, cmap=cmap)

# Cone of Interference のプロット
ax_sc.fill_between(tm, COI, fc='w', hatch='x', alpha=0.5)
ax_sc.plot(tm, COI, '--', color='black')

# 右上の文字
ax_sc.text(0.99, 0.97, "scalogram", color='white', ha='right', va='top',
              path_effects=[path_effects.Stroke(linewidth=2, foreground='black'),
                            path_effects.Normal()], 
              transform=ax_sc.transAxes)

mpl.colorbar.ColorbarBase(cx_sc, cmap=cmap,
                          norm=norm,
                          orientation="vertical",
                          label='log amplitude')

# 図を保存
plt.savefig(out_fig_path, transparent=False)

実際にコードを実行してみよう

上で記述したコードは cwt-demo.pyという名前のファイルとして,好きなディレクトリ(ここでは~/Desktop/LabCode/python/cwtとします)に保存します。
ターミナルを開き,

$ cd ~/Desktop/LabCode/python/cwt

と入力し,ディレクトリを移動します。あとは以下のコマンドで cwt-demo.pyを実行するだけです。( $マークは無視してください)。

$ python3 cwt-demo.py

実行結果

~/Desktop/LabCode/python/cwt にscalogram.pngというファイルができて,以下のようになっていれば成功です。

一段目は周波数と振幅が時間的に変化するテスト入力信号,二段目はスカログラムです。周波数と振幅の時間変化を追うことができていると思います!

右下と左下にみえるハッチをしてある部分はCOIの部分で,この部分はデータが有限であるためにスカログラムがエッジの部分の影響を受けていることを示しています。この部分の信頼度は低いです。

また,高周波の領域は低周波の領域と比べて上下に大きく広がっており,周波数分解能が悪いことがわかります(ただし時間分解能は良くなっています)。

コードの解説

作成したコードの解説をしていきます。

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patheffects as path_effects
import numpy as np
```

必要なモジュールをインポートします。

```python
#------------------------------------------------------------------------------
# テスト信号の生成
#------------------------------------------------------------------------------
dt = 1.0e-4
tms = 0.0
tme = 2.0
tm = np.arange(tms, tme, dt)

amp = 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))

テスト信号を作成します。振幅と振動数が時間的に変化するものを用います。ノイズも加えてより現実的にしておきます。

#------------------------------------------------------------------------------
# 諸々の設定
freq_upper = 300.0 # in kHz
out_fig_path = "scalogram.png" 

スカログラムの周波数上限の設定や出力される図の名前の設定です。

#------------------------------------------------------------------------------
# 連続ウェーブレット変換
# --- n_cwt をlen(sig)よりも大きな2のべき乗の数になるように設定
n_cwt = int(2**(np.ceil(np.log2(len(sig)))))

# --- 後で使うパラメータを定義
dj = 0.125
omega0 = 6.0
s0 = 2.0*dt
J = int(np.log2(n_cwt*dt/s0)/dj)

# --- スケール
s = s0*2.0**(dj*np.arange(0, J+1, 1))

# --- n_cwt個のデータになるようにゼロパディングをして,DC成分を除く
x = np.zeros(n_cwt)
x[0:len(sig)] = sig[0:len(sig)] - np.mean(sig)

# --- omega array
omega = 2.0*np.pi*np.fft.fftfreq(n_cwt, dt)

# --- FFTを使って離散ウェーブレット変換する
X = np.fft.fft(x)
cwt = np.zeros((J+1, n_cwt), dtype=complex) # CWT array 

Hev = np.array(omega > 0.0)
for j in range(J+1):
    Psi = np.sqrt(2.0*np.pi*s[j]/dt)*np.pi**(-0.25)*np.exp(-(s[j]*omega-omega0)**2/2.0)*Hev
    cwt[j, :] = np.fft.ifft(X*np.conjugate(Psi))

s_to_f = (omega0 + np.sqrt(2 + omega0**2)) / (4.0*np.pi) 
freq_cwt = s_to_f / s
cwt = cwt[:, 0:len(sig)]

# --- cone of interference
COI = np.zeros_like(tm)
COI[0] = 0.5/dt
COI[1:len(tm)//2] = np.sqrt(2)*s_to_f/tm[1:len(tm)//2]
COI[len(tm)//2:-1] = np.sqrt(2)*s_to_f/(tm[-1]-tm[len(tm)//2:-1])
COI[-1] = 0.5/dt

ウェーブレット変換の部分です。上で説明したように畳み込み定理を用いてFFTを使って計算しています。

#------------------------------------------------------------------------------
# 解析結果の可視化
figsize = (210/25.4, 294/25.4)
dpi = 200
fig = plt.figure(figsize=figsize, dpi=dpi)

# --- 図の設定 (全体)
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['xtick.top'] = True
plt.rcParams['xtick.major.size'] = 6
plt.rcParams['xtick.minor.size'] = 3
plt.rcParams['xtick.minor.visible'] = True
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['ytick.right'] = True
plt.rcParams['ytick.major.size'] = 6
plt.rcParams['ytick.minor.size'] = 3
plt.rcParams['ytick.minor.visible'] = True
plt.rcParams["font.size"] = 14
plt.rcParams['font.family'] = 'Helvetica'

# プロット枠 (axes) の設定
ax1 = fig.add_axes([0.15, 0.55, 0.70, 0.3])
ax_sc = fig.add_axes([0.15, 0.2, 0.70, 0.30])
cx_sc = fig.add_axes([0.87, 0.2, 0.02, 0.30])

# 元データのプロット
ax1.set_xlim(tms, tme)
ax1.set_xlabel('')
ax1.tick_params(labelbottom=False)
ax1.set_ylabel('x')

ax1.plot(tm, sig, c='black')

# スペクトログラムのプロット
ax_sc.set_xlim(tms, tme)
ax_sc.set_xlabel('time (s)')
ax_sc.tick_params(labelbottom=True)
ax_sc.set_ylim(0, freq_upper)
ax_sc.set_ylabel('frequency\n(Hz)')

# カラーバーのレンジの設定
norm = mpl.colors.Normalize(vmin=np.log10(np.abs(cwt[freq_cwt < freq_upper, :])**2).min(), 
                            vmax=np.log10(np.abs(cwt[freq_cwt < freq_upper, :])**2).max())
# カラーマップ
cmap = mpl.cm.jet

# スカログラムのプロット
ax_sc.contourf(tm, freq_cwt, np.log10(np.abs(cwt)**2), 
                norm=norm, levels=256, cmap=cmap)

# Cone of Interference のプロット
ax_sc.fill_between(tm, COI, fc='w', hatch='x', alpha=0.5)
ax_sc.plot(tm, COI, '--', color='black')

# 右上の文字
ax_sc.text(0.99, 0.97, "scalogram", color='white', ha='right', va='top',
              path_effects=[path_effects.Stroke(linewidth=2, foreground='black'),
                            path_effects.Normal()], 
              transform=ax_sc.transAxes)

mpl.colorbar.ColorbarBase(cx_sc, cmap=cmap,
                          norm=norm,
                          orientation="vertical",
                          label='log amplitude')

# 図を保存
plt.savefig(out_fig_path, transparent=False)

解析結果を可視化しています。この部分はSTFTのものとほとんど同じです。
STFTには無いCOIの部分は, fill_between を用いてCOI以下の部分のハッチを描画を,通常のplotを用いてCOI曲線の描画を行っています。

最後に

今回は,時間周波数解析の手法の一つであるウェーブレット解析について紹介しました。STFTと比べて窓関数幅の選択の必要がないなどの利点があります。
解析の目的,信号の特性に合わせて,適切な手法を選べるようにしましょう!

参考文献

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です