ウェーブレット変換入門：時間-周波数解析のPython実装

連続ウェーブレット変換（CWT）と離散ウェーブレット変換（DWT）の数学的基礎を解説し、PyWaveletsを使った時間-周波数解析のPython実装を紹介します。

4 min read

はじめに

高速フーリエ変換（FFT）の記事で解説したように、フーリエ変換は信号をさまざまな周波数の正弦波に分解する強力な手法です。しかし、フーリエ変換には本質的な制約があります。それは時間情報の喪失です。

フーリエ変換の結果は「信号全体にどの周波数成分がどれだけ含まれるか」を示しますが、「その周波数成分がいつ現れたか」はわかりません。音声信号や地震波のように周波数が時間とともに変化する信号に対して、この制約は致命的です。

短時間フーリエ変換（STFT）は窓関数を用いて時間局在性を導入しますが、窓サイズが固定であるため、高周波と低周波を同時に高い分解能で捉えることができません。ウェーブレット変換はこの問題を解決し、時間と周波数の両方において適応的な分解能を実現します。

本記事では、連続ウェーブレット変換（CWT）と離散ウェーブレット変換（DWT）の数学的基礎を解説し、PyWaveletsを用いたPython実装を紹介します。

フーリエ変換からウェーブレット変換へ

STFTの限界：固定窓サイズの問題

STFTは信号を短い窓で切り出し、各区間にフーリエ変換を適用することで時間-周波数表現を得ます。

\[\text{STFT}\{x(t)\}(\tau, \omega) = \int_{-\infty}^{\infty} x(t) \, w(t - \tau) \, e^{-j\omega t} \, dt \tag{1}\]

ここで \(w(t)\) は窓関数、\(\tau\) は時間シフトです。STFTの問題は、窓サイズが固定である点にあります。

  • 窓が短い場合：時間分解能は高いが、周波数分解能が低い
  • 窓が長い場合：周波数分解能は高いが、時間分解能が低い

このトレードオフはハイゼンベルクの不確定性原理に起因します。時間幅 \(\Delta t\) と周波数幅 \(\Delta f\) の間には次の下限が存在します。

\[\Delta t \cdot \Delta f \geq \frac{1}{4\pi} \tag{2}\]

STFTでは \(\Delta t\) と \(\Delta f\) の積が一定かつ固定ですが、ウェーブレット変換ではこの積を一定に保ちつつ、周波数に応じて \(\Delta t\) と \(\Delta f\) の配分を適応的に変化させます。

ウェーブレット変換のアイデア

ウェーブレット変換の核心は、正弦波の代わりに**局在化した波形（ウェーブレット）**を基底関数として用いることです。ウェーブレットは時間的に有限の広がりを持つため、時間局在性を自然に備えています。

さらに、ウェーブレットを**伸縮（スケーリング）**させることで、異なる周波数帯域を異なる分解能で解析できます。

  • 小さいスケール（圧縮されたウェーブレット）→ 高周波成分を高い時間分解能で捕捉
  • 大きいスケール（引き伸ばされたウェーブレット）→ 低周波成分を高い周波数分解能で捕捉

連続ウェーブレット変換（CWT）

マザーウェーブレット

ウェーブレット変換の基底となる関数をマザーウェーブレット \(\psi(t)\) と呼びます。マザーウェーブレットは次の2つの条件を満たす必要があります。

  1. 有限エネルギー：\(\int_{-\infty}^{\infty} |\psi(t)|^2 \, dt < \infty\)
  2. 許容条件（Admissibility condition）：\(C_\psi = \int_0^{\infty} \frac{|\hat{\Psi}(\omega)|^2}{\omega} \, d\omega < \infty\)

許容条件は、マザーウェーブレットの平均値がゼロであること（\(\int \psi(t) \, dt = 0\)）を意味します。つまり、ウェーブレットは振動的な波形でなければなりません。

CWTの定義

信号 \(x(t)\) の連続ウェーブレット変換は、スケールパラメータ \(a > 0\) と平行移動パラメータ \(b\) を用いて次のように定義されます。

\[W(a, b) = \frac{1}{\sqrt{a}} \int_{-\infty}^{\infty} x(t) \, \psi^*\!\left(\frac{t - b}{a}\right) dt \tag{3}\]

ここで、

  • \(a\): スケールパラメータ（伸縮）。\(a\) が大きいほどウェーブレットが引き伸ばされ、低周波成分を捕捉する
  • \(b\): 平行移動パラメータ。信号上の解析位置を指定する
  • \(\psi^*\): マザーウェーブレットの複素共役
  • \(\frac{1}{\sqrt{a}}\): エネルギー正規化係数。スケールが変わってもウェーブレットのエネルギーが一定に保たれる

CWTの結果 \(W(a, b)\) は2次元の時間-スケール表現であり、これをスカログラム（\(|W(a, b)|^2\)）として可視化します。

代表的なマザーウェーブレット

モルレーウェーブレット（Morlet wavelet）

\[\psi(t) = \pi^{-1/4} e^{j\omega_0 t} e^{-t^2/2} \tag{4}\]

ガウス窓で変調された複素正弦波であり、\(\omega_0\)（通常 \(\omega_0 = 5\) または \(6\)）は中心周波数を制御します。時間-周波数解析において最も広く使われるウェーブレットの一つで、周波数分解能と時間分解能のバランスに優れています。

メキシカンハットウェーブレット（Mexican Hat wavelet）

\[\psi(t) = \frac{2}{\sqrt{3}\pi^{1/4}} (1 - t^2) e^{-t^2/2} \tag{5}\]

ガウス関数の2階微分（負の符号付き）であり、ガウス差分（DoG: Difference of Gaussians）とも呼ばれます。実数値のウェーブレットで、信号の特異点やエッジの検出に適しています。

離散ウェーブレット変換（DWT）

ダイアディックサンプリング

CWTはスケールと平行移動を連続的に変化させるため、冗長な情報を多く含みます。計算効率と理論的な簡潔さを両立させるため、スケールと平行移動を2の冪乗で離散化したものが**離散ウェーブレット変換（DWT）**です。

\[a = 2^j, \quad b = k \cdot 2^j \tag{6}\]

ここで \(j\) は分解レベル（スケールレベル）、\(k\) は平行移動インデックスです。これにより、ウェーブレット基底は次のようになります。

\[\psi_{j,k}(t) = 2^{-j/2} \psi(2^{-j} t - k) \tag{7}\]

多重解像度解析（MRA）

DWTの理論的基盤となるのが多重解像度解析（MRA: Multiresolution Analysis）です。MRAでは、信号を異なる解像度レベルの近似成分と詳細成分に分解します。

MRAは2つの関数によって構成されます。

  • スケーリング関数 \(\phi(t)\)：低周波成分（近似）を表現する
  • ウェーブレット関数 \(\psi(t)\)：高周波成分（詳細）を表現する

スケーリング関数は次の再帰的関係（二重スケール方程式）を満たします。

\[\phi(t) = \sqrt{2} \sum_k h[k] \, \phi(2t - k) \tag{8}\]

ウェーブレット関数はスケーリング関数から導かれます。

\[\psi(t) = \sqrt{2} \sum_k g[k] \, \phi(2t - k) \tag{9}\]

ここで \(h[k]\) はローパスフィルタ係数、\(g[k]\) はハイパスフィルタ係数であり、次の関係があります。

\[g[k] = (-1)^k h[N - 1 - k] \tag{10}\]

フィルタバンクによる分解

DWTの計算は、フィルタバンクとダウンサンプリングの繰り返しとして効率的に実装できます。これはMallatのアルゴリズムとして知られています。

1段の分解では、信号をローパスフィルタ \(H\) とハイパスフィルタ \(G\) に通し、それぞれの出力を2分の1にダウンサンプリングします。

  • 近似係数（cA: Approximation Coefficients）：ローパスフィルタ \(H\) + ダウンサンプリング → 低周波成分
  • 詳細係数（cD: Detail Coefficients）：ハイパスフィルタ \(G\) + ダウンサンプリング → 高周波成分

多レベル分解では、近似係数 cA をさらに再帰的に分解していきます。レベル \(j\) の分解では、信号は次の成分に分離されます。

\[x(t) = \text{cA}_j + \text{cD}_j + \text{cD}_{j-1} + \cdots + \text{cD}_1 \tag{11}\]

この構造は、ローパスフィルタの設計と比較で扱ったフィルタの概念をカスケード接続したものと捉えることができます。

PyWaveletsによるPython実装

CWTの例：チャープ信号のスカログラム

周波数が時間とともに変化するチャープ信号に対してCWTを適用し、スカログラムを可視化します。

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

# --- チャープ信号の生成 ---
fs = 1000  # サンプリング周波数 [Hz]
t = np.arange(0, 2, 1 / fs)
# 周波数が10Hzから100Hzに線形に変化するチャープ信号
freq = np.linspace(10, 100, len(t))
signal = np.sin(2 * np.pi * np.cumsum(freq) / fs)

# --- CWTの計算 ---
scales = np.arange(1, 128)
coefficients, frequencies = pywt.cwt(signal, scales, "morl", 1 / fs)

# --- スカログラムの可視化 ---
plt.figure(figsize=(10, 5))
plt.pcolormesh(
    t, frequencies, np.abs(coefficients) ** 2, shading="gouraud", cmap="viridis"
)
plt.colorbar(label="Power")
plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")
plt.ylim(0, 150)
plt.title("CWT Scalogram (Morlet Wavelet)")
plt.tight_layout()
plt.show()

スカログラムでは、周波数が時間とともに10Hzから100Hzへ上昇していく様子が明確に可視化されます。STFTのスペクトログラムと比較すると、低周波域では周波数分解能が高く、高周波域では時間分解能が高いという適応的な分解能が確認できます。

DWTの例：多レベル分解

ノイズが重畳した信号に対してDWTによる多レベル分解を行い、各レベルの近似係数と詳細係数を可視化します。

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

# --- テスト信号の生成 ---
np.random.seed(42)
fs = 500
t = np.arange(0, 1, 1 / fs)
# 低周波成分 + 高周波成分 + ノイズ
clean = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 50 * t)
noise = 0.3 * np.random.randn(len(t))
signal = clean + noise

# --- DWTによる多レベル分解 ---
wavelet = "db4"  # Daubechies-4ウェーブレット
level = 5
coeffs = pywt.wavedec(signal, wavelet, level=level)
# coeffs = [cA5, cD5, cD4, cD3, cD2, cD1]

# --- 各レベルの可視化 ---
fig, axes = plt.subplots(level + 2, 1, figsize=(10, 10), sharex=False)

axes[0].plot(t, signal, "k", linewidth=0.5)
axes[0].set_title("Original Signal")
axes[0].set_ylabel("Amplitude")

axes[1].plot(coeffs[0], "b", linewidth=0.5)
axes[1].set_title(f"Approximation Coefficients (cA{level})")
axes[1].set_ylabel("Amplitude")

for i in range(1, level + 1):
    axes[i + 1].plot(coeffs[i], "r", linewidth=0.5)
    axes[i + 1].set_title(f"Detail Coefficients (cD{level - i + 1})")
    axes[i + 1].set_ylabel("Amplitude")

plt.tight_layout()
plt.show()

分解結果を観察すると、以下のことがわかります。

  • cA5（近似係数）：信号の最も低周波な成分を含む（5Hzの正弦波）
  • cD1, cD2（詳細係数、高レベル）：高周波ノイズが主に含まれる
  • cD3, cD4（詳細係数、中レベル）：50Hzの正弦波成分が含まれる

ウェーブレット閾値処理によるノイズ除去

DWTの重要な応用の一つが**ウェーブレット閾値処理（wavelet thresholding）**によるノイズ除去です。ノイズ成分は主に詳細係数に現れるため、閾値以下の係数をゼロにすることでノイズを除去できます。

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

# --- テスト信号の生成 ---
np.random.seed(42)
fs = 500
t = np.arange(0, 1, 1 / fs)
clean = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 50 * t)
noise = 0.3 * np.random.randn(len(t))
signal = clean + noise

# --- DWT分解 ---
wavelet = "db4"
level = 5
coeffs = pywt.wavedec(signal, wavelet, level=level)

# --- 閾値処理 ---
# ユニバーサル閾値（VisuShrink）: λ = σ * √(2 * log(N))
sigma = np.median(np.abs(coeffs[-1])) / 0.6745  # ノイズのロバスト推定
threshold = sigma * np.sqrt(2 * np.log(len(signal)))

# ソフト閾値処理を詳細係数に適用
coeffs_thresh = [coeffs[0]]  # 近似係数はそのまま
for c in coeffs[1:]:
    coeffs_thresh.append(pywt.threshold(c, threshold, mode="soft"))

# --- 逆DWTによる信号復元 ---
denoised = pywt.waverec(coeffs_thresh, wavelet)

# --- 結果の可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)

axes[0].plot(t, clean, "b", linewidth=0.8)
axes[0].set_title("Original Clean Signal")
axes[0].set_ylabel("Amplitude")

axes[1].plot(t, signal, "k", linewidth=0.5)
axes[1].set_title(f"Noisy Signal (SNR = {10*np.log10(np.var(clean)/np.var(noise)):.1f} dB)")
axes[1].set_ylabel("Amplitude")

axes[2].plot(t, denoised[: len(t)], "r", linewidth=0.8)
axes[2].set_title(f"Denoised Signal (Wavelet Thresholding)")
axes[2].set_ylabel("Amplitude")

for ax in axes:
    ax.grid(True, alpha=0.3)
axes[2].set_xlabel("Time [s]")

plt.tight_layout()
plt.show()

# --- 定量評価 ---
mse_noisy = np.mean((signal - clean) ** 2)
mse_denoised = np.mean((denoised[: len(t)] - clean) ** 2)
print(f"ノイズあり信号のMSE:   {mse_noisy:.6f}")
print(f"ノイズ除去後のMSE:     {mse_denoised:.6f}")
print(f"MSE改善率:             {(1 - mse_denoised / mse_noisy) * 100:.1f}%")

ソフト閾値処理では、閾値 \(\lambda\) を超える係数を \(\lambda\) だけ縮小し、閾値以下の係数をゼロにします。

\[\text{soft}(x, \lambda) = \text{sign}(x) \cdot \max(|x| - \lambda, 0) \tag{12}\]

ユニバーサル閾値 \(\lambda = \sigma\sqrt{2 \log N}\)（VisuShrink）は、ノイズの標準偏差 \(\sigma\) と信号長 \(N\) から自動的に決定されるため、パラメータ調整が不要です。\(\sigma\) はMAD（Median Absolute Deviation）を用いてロバストに推定しています。

FFT vs STFT vs ウェーブレット変換の比較

特性FFTSTFTウェーブレット変換
時間分解能なし固定適応的（高周波で高い）
周波数分解能高い（全体）固定適応的（低周波で高い）
基底関数正弦波窓付き正弦波マザーウェーブレット
出力1次元スペクトル2次元（時間-周波数）2次元（時間-スケール）
逆変換ありありあり
計算量\(O(N \log N)\)\(O(N L \log L)\)CWT: \(O(N S)\), DWT: \(O(N)\)
主な用途定常信号の解析音声・音楽解析非定常信号・ノイズ除去

ここで \(L\) はSTFTの窓長、\(S\) はCWTのスケール数です。DWTはフィルタバンクによる実装で \(O(N)\) の計算量を達成しており、非常に効率的です。

各手法の使い分けの指針は次のとおりです。

  • 信号が定常的で全体の周波数成分を知りたい → FFT
  • 時間変化する周波数を一様な分解能で追跡したい → STFT
  • 低周波の微細な変化と高周波の瞬時的なイベントを同時に捉えたい → ウェーブレット変換

まとめ

本記事では、ウェーブレット変換の数学的基礎からPython実装までを解説しました。

  • CWTはスケールと平行移動を連続的に変化させ、適応的な時間-周波数分解能を実現する
  • DWTはダイアディックサンプリングとフィルタバンクにより \(O(N)\) の効率的な分解を実現する
  • ウェーブレット閾値処理は、詳細係数に閾値を適用することでノイズ除去を行う手法として実用的
  • FFT、STFT、ウェーブレット変換はそれぞれ得意分野が異なり、解析対象に応じた使い分けが重要

ウェーブレット変換は、信号処理にとどまらず、画像圧縮（JPEG 2000）、金融データ解析、医療信号処理、地震波解析など、幅広い分野で活用されています。

関連記事

参考文献

関連記事