ウェーブレット変換の基礎理論と実装

フーリエ変換は信号に「どの周波数が含まれるか」を教えてくれますが、「その周波数がいつ出現するか」は教えてくれません。音楽を聴いているとき、「10秒目にドの音が鳴り、20秒目にレの音に変わった」という時間情報は、フーリエ変換のスペクトルからは失われてしまいます。

この時間-周波数の両方の情報を同時に捉えるのがウェーブレット変換（wavelet transform）です。短い「波のかけら」（wavelet）をスケーリングとシフトさせながら信号との類似度を測ることで、「いつ」「どのスケールの」成分が存在するかを明らかにします。

ウェーブレット変換を理解すると、以下のような応用が開けます。

信号処理: 非定常信号の時間-周波数解析
画像圧縮: JPEG 2000のコア技術
地震学: 地震波の到達時間と周波数の解析
金融工学: 株価の多重スケール解析
医療: 心電図（ECG）の異常検出

本記事の内容

フーリエ変換の限界と時間-周波数解析の必要性
連続ウェーブレット変換（CWT）の定義
マザーウェーブレットの種類（Morlet, Mexican hat等）
離散ウェーブレット変換（DWT）と多重解像度解析
Pythonによる実装と可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

フーリエ変換の限界

時間情報の喪失

フーリエ変換 $\hat{f}(\omega) = \int f(t) e^{-i\omega t}\, dt$ は信号全体にわたって積分するため、周波数情報は得られますが時間情報は失われます。

たとえば「前半に10Hz、後半に50Hz」の信号と「最初から最後まで10Hzと50Hzの混合」は、フーリエ変換では同じスペクトルになってしまいます。

身近な例で考えてみましょう。ピアノの演奏を録音した場合、「ド→レ→ミ」と弾いたのか「ド・レ・ミを同時に弾いた和音」なのかは、聴けば明らかに区別できます。しかしフーリエ変換のスペクトルを見ると、どちらの場合も同じ3つの周波数にピークが現れ、区別がつきません。楽譜に書かれた情報 — どの音がいつ鳴るか — がフーリエ変換では失われてしまうのです。

数学的に見ると、フーリエ変換の基底関数 $e^{-i\omega t}$ は $t = -\infty$ から $t = +\infty$ まで一様に広がる永遠の正弦波です。この基底関数と信号の内積を取ることは、信号の「全時刻にわたる平均的な周波数成分」を抽出する操作にほかなりません。時間に局在した情報は、この全時間平均の中に埋もれてしまうのです。

短時間フーリエ変換（STFT）の限界

時間情報を保持する素朴な方法として、信号を短い窓で切り出してからフーリエ変換する短時間フーリエ変換（STFT）があります。

$$ \text{STFT}(t, \omega) = \int f(\tau) w(\tau – t) e^{-i\omega\tau}\, d\tau $$

しかし窓関数 $w$ のサイズは固定なので、不確定性原理 $\Delta t \cdot \Delta\omega \geq 1/2$ のトレードオフが一律に適用されます。低周波を捉えるには長い窓が必要ですが、そうすると時間分解能が悪化します。

不確定性原理 — 時間と周波数の根本的なトレードオフ

不確定性原理（ガボール限界、ハイゼンベルクの不確定性原理のフーリエ解析版）は、任意の関数 $f(t)$ に対して

$$ \Delta t \cdot \Delta \omega \geq \frac{1}{2} $$

が成り立つことを述べています。ここで $\Delta t$ は信号の時間幅（時間方向の標準偏差）、$\Delta\omega$ は周波数幅（周波数方向の標準偏差）です。

これは「時間と周波数の両方を同時に無限に鋭く知ることはできない」という基本的な制約です。STFTでは窓幅を固定するため、$\Delta t$ と $\Delta\omega$ の比率が全周波数で一定になります。短い窓（小さな $\Delta t$）を使えば時間分解能は良くなりますが周波数分解能 $\Delta\omega$ が悪化し、長い窓を使えばその逆になります。

具体的な数値で考えましょう。サンプリング周波数1,000Hzの信号をSTFTで解析する場合、窓幅を100ms（100点）に設定すると、周波数分解能は $\Delta f = 1/0.1 = 10$ Hzです。5Hzと15Hzの正弦波は区別できますが、50Hzと55Hzの短いバースト信号が何ミリ秒に出現したかは、100msの窓の中では曖昧になります。窓幅を10msに短くすれば時間分解能は改善しますが、周波数分解能が $\Delta f = 100$ Hzに悪化して、低周波成分を全く分解できなくなります。

ウェーブレット変換はこの問題を解決します。スケール（周波数に対応）に応じて窓の幅を自動的に変えるのです。低周波には広い窓を使って周波数分解能を確保し、高周波には狭い窓を使って時間分解能を確保します。数学的には、不確定性原理の下限 $\Delta t \cdot \Delta \omega = 1/2$ に近い最適な配分を各スケールで実現するのがウェーブレット変換の本質です。

フーリエ変換とSTFTの限界を踏まえたうえで、これらの問題を解決するウェーブレット変換の定義を見ていきましょう。

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

定義

信号 $f(t)$ の連続ウェーブレット変換は

$$ \begin{equation} W_f(a, b) = \frac{1}{\sqrt{|a|}}\int_{-\infty}^{\infty} f(t)\, \psi^*\left(\frac{t – b}{a}\right) dt \end{equation} $$

ここで

$\psi(t)$: マザーウェーブレット（母ウェーブレット）— 局在した波形
$a$: スケールパラメータ — $a$ が大きいと波が広がり（低周波）、小さいと縮む（高周波）
$b$: シフトパラメータ — 時間方向の位置
$1/\sqrt{|a|}$: エネルギー正規化

$W_f(a, b)$ は「スケール $a$、位置 $b$ での信号とウェーブレットの類似度」を表します。

CWTの直感的理解 — 音叉を信号に当てるアナロジー

CWTの操作を直感的に理解するために、音叉（おんさ）を使ったアナロジーを考えましょう。

音叉はそれぞれ固有の周波数で振動します。440Hzの音叉を鳴っている楽器に近づけると、440Hz成分が含まれていれば共鳴して振動が大きくなります。CWTとは、「さまざまな周波数の音叉を、信号の各時刻に順番に当てていく操作」と考えることができます。

スケール $a$ を変えることは、音叉の固有周波数を変えることに対応します。$a$ を大きくするとウェーブレットが引き伸ばされて低周波の音叉になり、$a$ を小さくすると高周波の音叉になります
シフト $b$ を変えることは、音叉を信号上の異なる時刻に当てることに対応します
$W_f(a, b)$ の大きさは、その時刻・その周波数で信号がどれだけ「共鳴」したかを表します

スケールと周波数の関係

スケールパラメータ $a$ と対応する擬似周波数 $f_a$ の関係は

$$ f_a = \frac{f_c}{a \cdot \Delta t} $$

で与えられます。ここで $f_c$ はマザーウェーブレットの中心周波数、$\Delta t$ はサンプリング間隔です。Morletウェーブレット（$\omega_0 = 5$）の場合 $f_c = \omega_0 / (2\pi) \approx 0.796$ です。

この逆比例関係により、$a = 1$（小スケール）は高周波に、$a = 100$（大スケール）は低周波に対応します。

エネルギー正規化の意味

定義式の前にある $1/\sqrt{|a|}$ の正規化因子は、スケールを変えてもウェーブレットのエネルギー（$L^2$ ノルム）を一定に保つために必要です。スケール $a$ で引き伸ばされたウェーブレット $\psi(t/a)$ のエネルギーは

$$ \int |\psi(t/a)|^2\, dt = |a| \int |\psi(s)|^2\, ds = |a| \|\psi\|^2 $$

$s = t/a$ の変数変換で $dt = |a| ds$ となるためです。$1/\sqrt{|a|}$ を掛けることでこのスケール依存性を打ち消し、大スケールと小スケールのウェーブレット係数を公平に比較できるようにします。

マザーウェーブレットの条件

$\psi(t)$ がウェーブレットとして使えるための条件（許容条件）は

$$ C_\psi = \int_0^{\infty} \frac{|\hat{\psi}(\omega)|^2}{\omega}\, d\omega < \infty $$

この条件は $\hat{\psi}(0) = 0$、すなわち $\int \psi(t)\, dt = 0$（ゼロ平均）を必要とします。つまりウェーブレットは上下に振動する波形でなければなりません。

許容条件が成り立つ理由を直感的に説明すると、$C_\psi$ が有限であることは、周波数ゼロ（$\omega = 0$）付近でのスペクトルのエネルギーが急速にゼロに近づくことを要求しています。$|\hat{\psi}(\omega)|^2 / \omega$ が $\omega \to 0$ で発散しないためには、$\hat{\psi}(0) = 0$ が必要です。これは $\psi$ がゼロ平均、すなわち正の部分と負の部分が完全に打ち消し合うことを意味します。

ガウス関数 $e^{-t^2/2}$ はゼロ平均ではないため（常に正）、そのままではウェーブレットとして使えません。しかし、ガウス関数の2階微分（Mexican hat）やガウス関数に正弦波を掛けたもの（Morlet）はゼロ平均になるため、ウェーブレットの条件を満たします。

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

Morletウェーブレット:

$$ \psi(t) = e^{i\omega_0 t} e^{-t^2/2} $$

複素正弦波にガウス窓を掛けたもので、$\omega_0 = 5$ 程度が標準です。周波数分析に最も広く使われます。Morletウェーブレットが優れている点は、時間-周波数平面でのガウス型の局在性にあります。不確定性原理の下限に最も近い（$\Delta t \cdot \Delta \omega \approx 1/2$）性質を持ち、時間分解能と周波数分解能のバランスが最適に近いです。

Morletウェーブレットのフーリエ変換は

$$ \hat{\psi}(\omega) = \pi^{-1/4} e^{-(\omega – \omega_0)^2/2} $$

となり、$\omega = \omega_0$ を中心とするガウス関数です。したがってスケール $a$ でのフーリエ変換は $\omega = \omega_0/a$ にピークを持ち、対応する周波数が $f = \omega_0 / (2\pi a)$ であることが確認できます。

Mexican hatウェーブレット（Ricker wavelet）:

$$ \psi(t) = \frac{2}{\sqrt{3}\pi^{1/4}}(1 – t^2) e^{-t^2/2} $$

ガウス関数の負の正規化2階微分に比例します。その名の通り、形がメキシコの帽子に似ています。実数値で、エッジや特異点の検出に適しています。地震学では地震波形の解析にRicker waveletが標準的に用いられています。

Mexican hatのフーリエ変換は

$$ \hat{\psi}(\omega) = \frac{2}{\sqrt{3}\pi^{1/4}} \sqrt{2\pi}\, \omega^2 e^{-\omega^2/2} $$

であり、$\omega = \sqrt{2}$ にピークを持ちます。

Haarウェーブレット:

$$ \psi(t) = \begin{cases} 1 & (0 \leq t < 1/2) \\ -1 & (1/2 \leq t < 1) \\ 0 & (\text{otherwise}) \end{cases} $$

最も単純なウェーブレットであり、1909年にAlfred Haarによって導入された歴史的にも最初のウェーブレットです。不連続であるため周波数分解能は低いですが、計算が極めて高速であり、離散ウェーブレット変換の基本として重要な役割を果たします。Haarウェーブレットの係数は「隣接する2点の差」に対応するため、信号の急激な変化（ステップ変化）を検出するのに適しています。

ウェーブレットの選択は解析の目的に依存します。一般的な指針として、周波数分析にはMorlet、エッジ検出にはMexican hat、高速計算や画像圧縮にはDaubechiesウェーブレット（Haarの一般化）が使われます。

各マザーウェーブレットの特徴を理解したところで、CWTの逆変換と復元可能性について見ていきましょう。

逆変換

$$ f(t) = \frac{1}{C_\psi}\int_0^{\infty}\int_{-\infty}^{\infty} W_f(a,b)\, \frac{1}{\sqrt{a}}\psi\left(\frac{t-b}{a}\right) \frac{db\, da}{a^2} $$

CWTは可逆であり、ウェーブレット係数 $W_f(a, b)$ から元の信号 $f(t)$ を完全に復元できます。$C_\psi$ は許容条件で定義された定数であり、この定数が有限であること（許容条件が成り立つこと）が逆変換の存在を保証します。

逆変換の直感的な意味は、全てのスケール $a$ と位置 $b$ のウェーブレット $\psi_{a,b}(t)$ に対応する係数 $W_f(a,b)$ を重み付きで足し合わせることで、元の信号を「再合成」する操作です。フーリエ逆変換が各周波数成分を足し合わせて元の信号を復元するのと同様に、CWTの逆変換は各スケール・各位置のウェーブレット成分を足し合わせます。

積分の測度 $db\, da / a^2$ は、スケール軸上での「面積要素」であり、小さなスケール（高周波）にも大きなスケール（低周波）にも等しい重みを与えるためのものです。

連続ウェーブレット変換は理論的に美しいですが、計算量が膨大（全てのスケールとシフトを連続的に計算）であるため、実用的には離散化が不可欠です。次に、効率的な計算を可能にする離散ウェーブレット変換を見ていきましょう。

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

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

CWTは $(a, b)$ が連続なので冗長です。実用的には $a = 2^j$, $b = k \cdot 2^j$ とダイアディック（2のべき乗）にサンプリングします。

$$ \psi_{j,k}(t) = 2^{-j/2}\psi(2^{-j}t – k) $$

ここで $j$ はスケールレベル（$j$ が大きいほど低周波）、$k$ はシフト位置を表す整数です。$2^{-j/2}$ はエネルギー正規化因子です。

なぜ2のべき乗でサンプリングするのでしょうか。これは「オクターブ」の概念に対応します。音楽において1オクターブ上の音は周波数が2倍です。同様に、各スケールレベルで周波数帯域が2倍ずつ変化するダイアディックサンプリングは、対数的に等間隔な周波数分割を実現します。この分割は人間の聴覚特性（メル尺度）にも近く、知覚的にも自然な分解方法です。

ダイアディックサンプリングの重要な性質は、適切なウェーブレットを選べば $\{\psi_{j,k}\}$ が $L^2(\mathbb{R})$ の正規直交基底を構成することです。これにより、情報の損失なしに効率的な係数表現が可能になります。

フィルタバンク解釈

DWTはローパスフィルタとハイパスフィルタの組（フィルタバンク）として効率的に実装できます。この等価性はStephane Mallat（マラー）が1989年に示したもので、Mallatのアルゴリズムと呼ばれています。

信号 $x[n]$ に対して

ローパスフィルタ $h[n]$ で畳み込み → ダウンサンプリング（2点おきに取る） → 近似係数 $a[n]$
ハイパスフィルタ $g[n]$ で畳み込み → ダウンサンプリング → 詳細係数 $d[n]$

近似係数 $a[n]$ は信号の「大まかな形状」（低周波成分）を表し、詳細係数 $d[n]$ は「細かい変動」（高周波成分）を表します。

この操作を近似係数に対して再帰的に繰り返すのが多重解像度解析（Multi-Resolution Analysis, MRA）です。

Haarウェーブレットでの具体的な計算例

Haarウェーブレットの場合、ローパスフィルタは $h = [1/\sqrt{2},\, 1/\sqrt{2}]$、ハイパスフィルタは $g = [1/\sqrt{2},\, -1/\sqrt{2}]$ です。

入力信号 $x = [4, 6, 10, 2, 8, 4, 6, 2]$（長さ8）に対して1レベルのDWTを計算しましょう。

近似係数（隣接2点の平均の $\sqrt{2}$ 倍）:

$$ \begin{aligned} a[0] &= (x[0] + x[1])/\sqrt{2} = (4 + 6)/\sqrt{2} = 10/\sqrt{2} \approx 7.07 \\ a[1] &= (x[2] + x[3])/\sqrt{2} = (10 + 2)/\sqrt{2} = 12/\sqrt{2} \approx 8.49 \\ a[2] &= (x[4] + x[5])/\sqrt{2} = (8 + 4)/\sqrt{2} = 12/\sqrt{2} \approx 8.49 \\ a[3] &= (x[6] + x[7])/\sqrt{2} = (6 + 2)/\sqrt{2} = 8/\sqrt{2} \approx 5.66 \end{aligned} $$

詳細係数（隣接2点の差の $\sqrt{2}$ 分の1倍）:

$$ \begin{aligned} d[0] &= (x[0] – x[1])/\sqrt{2} = (4 – 6)/\sqrt{2} = -2/\sqrt{2} \approx -1.41 \\ d[1] &= (x[2] – x[3])/\sqrt{2} = (10 – 2)/\sqrt{2} = 8/\sqrt{2} \approx 5.66 \\ d[2] &= (x[4] – x[5])/\sqrt{2} = (8 – 4)/\sqrt{2} = 4/\sqrt{2} \approx 2.83 \\ d[3] &= (x[6] – x[7])/\sqrt{2} = (6 – 2)/\sqrt{2} = 4/\sqrt{2} \approx 2.83 \end{aligned} $$

元の信号（8点）が近似係数（4点）+ 詳細係数（4点）= 8つの係数に変換されています。データ量は変わりませんが、信号の低周波成分と高周波成分が分離されています。$d[1] \approx 5.66$ が最も大きいのは、$x[2] = 10$ と $x[3] = 2$ の間に大きな変動があることを反映しています。

完全再構成条件

DWTが有用であるためには、分解した係数から元の信号を完全に復元できる必要があります。これを保証するのが完全再構成条件です。フィルタ $h$ と $g$ が以下を満たすとき、完全再構成が可能です。

$$ |H(\omega)|^2 + |G(\omega)|^2 = 2 \quad (\text{for all } \omega) $$

ここで $H(\omega)$ と $G(\omega)$ はそれぞれ $h[n]$ と $g[n]$ のフーリエ変換です。Haarフィルタの場合、$H(\omega) = (1 + e^{-i\omega})/\sqrt{2}$ なので $|H(\omega)|^2 = (1 + \cos\omega)$ であり、$|G(\omega)|^2 = (1 – \cos\omega)$ と合わせて和が2になることが確認できます。

計算量の優位性

DWTの計算量は $O(N)$ です。これはFFTの $O(N \log N)$ よりも高速です。各レベルでの計算量が $O(N_j)$ であり、$N_j = N/2^j$ なので総計算量は

$$ \sum_{j=1}^{J} O(N/2^j) = O(N) \sum_{j=1}^{J} 2^{-j} < O(N) \cdot 1 = O(N) $$

等比級数の和が2未満に収束するためです。この計算効率の高さがDWTの実用上の大きな利点です。

DWTの理論的な枠組みを理解したところで、PythonでCWTとDWTを実装して実際の信号を解析してみましょう。

Pythonでの実装

CWTの実装と可視化

import numpy as np
import matplotlib.pyplot as plt

def morlet_wavelet(t, omega0=5.0):
    """Morletウェーブレット"""
    return np.exp(1j * omega0 * t) * np.exp(-t**2 / 2) * np.pi**(-0.25)

def cwt(signal, scales, dt=1.0, omega0=5.0):
    """連続ウェーブレット変換"""
    N = len(signal)
    t = np.arange(N) * dt
    coeffs = np.zeros((len(scales), N), dtype=complex)

    for i, a in enumerate(scales):
        # スケーリングされたウェーブレット
        t_scaled = (t - t[N//2]) / a
        wavelet = morlet_wavelet(t_scaled, omega0) / np.sqrt(a)
        # 畳み込み（FFTで高速化）
        coeffs[i] = np.convolve(signal, np.conj(wavelet[::-1]), mode="same") * dt

    return coeffs

# テスト信号: 時間とともに周波数が変化するチャープ信号
np.random.seed(42)
fs = 500
t = np.arange(0, 4, 1/fs)
N = len(t)

# 前半: 10Hz、中間: 30Hz、後半: 10Hz + 50Hz
signal = np.zeros(N)
signal[:N//3] = np.sin(2*np.pi*10*t[:N//3])
signal[N//3:2*N//3] = np.sin(2*np.pi*30*t[N//3:2*N//3])
signal[2*N//3:] = np.sin(2*np.pi*10*t[2*N//3:]) + np.sin(2*np.pi*50*t[2*N//3:])
signal += 0.2 * np.random.randn(N)

# CWT計算
scales = np.arange(1, 100)
freqs = 5.0 * fs / (2 * np.pi * scales)  # Morletの周波数変換
coeffs = cwt(signal, scales, dt=1/fs)

fig, axes = plt.subplots(3, 1, figsize=(14, 10))

# (a) 時間波形
ax = axes[0]
ax.plot(t, signal, "b-", linewidth=0.5)
ax.set_ylabel("Amplitude", fontsize=12)
ax.set_title("Signal: Frequency changes over time", fontsize=13)
ax.grid(True, alpha=0.3)

# (b) スカログラム（CWTの絶対値）
ax = axes[1]
im = ax.pcolormesh(t, freqs, np.abs(coeffs), shading="auto", cmap="hot")
ax.set_ylabel("Frequency (Hz)", fontsize=12)
ax.set_ylim(0, 80)
ax.set_title("CWT Scalogram (Morlet wavelet)", fontsize=13)
plt.colorbar(im, ax=ax, label="|W(a,b)|")

# (c) STFTとの比較（スペクトログラム）
ax = axes[2]
from scipy.signal import spectrogram
f_stft, t_stft, Sxx = spectrogram(signal, fs, nperseg=256, noverlap=200)
ax.pcolormesh(t_stft, f_stft, 10*np.log10(Sxx+1e-10), shading="auto", cmap="hot")
ax.set_ylabel("Frequency (Hz)", fontsize=12)
ax.set_xlabel("Time (s)", fontsize=12)
ax.set_ylim(0, 80)
ax.set_title("STFT Spectrogram (fixed window)", fontsize=13)

plt.tight_layout()
plt.savefig("wavelet_cwt.png", dpi=150, bbox_inches="tight")
plt.show()

このグラフから、ウェーブレット変換の時間-周波数解析能力が確認できます。

上段（信号）: 前半は10Hz、中盤は30Hz、後半は10Hzと50Hzの合成という、時間とともに周波数構成が変化する非定常信号です。
中段（CWTスカログラム）: Morletウェーブレットによるスカログラムが、各時間区間での周波数成分を正確に捉えています。低周波（10Hz）では時間分解能がやや粗いですが周波数分解能が高く、高周波（50Hz）では時間分解能が良いことが確認できます。これがウェーブレット変換の「適応的な分解能」の本質です。
下段（STFTスペクトログラム）: 固定窓のSTFTでも時間-周波数の変化を捉えていますが、窓サイズが固定であるため、低周波と高周波で同じ分解能になっています。ウェーブレット変換の方がスケールに応じた適応的な分解能を持つ点で優れています。

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

import numpy as np
import matplotlib.pyplot as plt

def haar_dwt_1level(x):
    """Haarウェーブレットによる1レベルDWT"""
    N = len(x)
    approx = (x[0::2] + x[1::2]) / np.sqrt(2)  # 近似係数
    detail = (x[0::2] - x[1::2]) / np.sqrt(2)  # 詳細係数
    return approx, detail

def haar_dwt_multilevel(x, levels):
    """多重レベルDWT"""
    details = []
    current = x.copy()
    for _ in range(levels):
        approx, detail = haar_dwt_1level(current)
        details.append(detail)
        current = approx
    return current, details

# テスト信号
np.random.seed(42)
N = 1024
t = np.linspace(0, 1, N)
# 段階的に周波数が変わる信号
signal = np.zeros(N)
signal[:256] = np.sin(2*np.pi*5*t[:256])
signal[256:512] = np.sin(2*np.pi*20*t[256:512])
signal[512:768] = np.sin(2*np.pi*50*t[512:768])
signal[768:] = np.sin(2*np.pi*100*t[768:])

# 4レベルDWT
levels = 4
approx, details = haar_dwt_multilevel(signal, levels)

fig, axes = plt.subplots(levels + 2, 1, figsize=(14, 12), sharex=False)

# 元の信号
ax = axes[0]
ax.plot(t, signal, "b-", linewidth=0.5)
ax.set_title("Original Signal", fontsize=12)
ax.set_ylabel("Amp", fontsize=10)
ax.grid(True, alpha=0.3)

# 詳細係数（各レベル）
for i, d in enumerate(details):
    ax = axes[i + 1]
    t_d = np.linspace(0, 1, len(d))
    ax.plot(t_d, d, "r-", linewidth=0.5)
    ax.set_title(f"Detail coefficients (Level {i+1}, scale $2^{{{i+1}}}$)", fontsize=11)
    ax.set_ylabel("Amp", fontsize=10)
    ax.grid(True, alpha=0.3)

# 近似係数（最終レベル）
ax = axes[-1]
t_a = np.linspace(0, 1, len(approx))
ax.plot(t_a, approx, "g-", linewidth=1)
ax.set_title(f"Approximation coefficients (Level {levels})", fontsize=11)
ax.set_xlabel("Time (normalized)", fontsize=12)
ax.set_ylabel("Amp", fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("wavelet_dwt.png", dpi=150, bbox_inches="tight")
plt.show()

このグラフから、DWTによる多重解像度解析の仕組みを詳しく読み取ることができます。

元の信号は4つの区間でそれぞれ異なる周波数（5Hz, 20Hz, 50Hz, 100Hz）を持つ非定常信号です。DWTの各レベルの詳細係数は異なる周波数帯域に対応しています。

Level 1（最も高い周波数帯域）: 100Hzの信号が存在する $t = 0.75$〜$1.0$ の区間で詳細係数の振幅が大きくなっています。サンプリング周波数を $f_s$ とすると、Level 1 は $f_s/4$ 〜 $f_s/2$ の帯域を担当します
Level 2: 50Hzの信号が存在する $t = 0.5$〜$0.75$ の区間で振幅が大きくなっています。$f_s/8$ 〜 $f_s/4$ の帯域に対応します
Level 3: 20Hzの信号がある $t = 0.25$〜$0.5$ 区間で反応しています
Level 4: 最も低い周波数帯域を担当し、5Hzの成分を捉えています

最下段の近似係数（Approximation）は、Level 4 の分解後に残った最低周波数成分であり、信号の大域的なトレンドを表しています。

注目すべきは、各レベルの係数のデータ点数が半分ずつ減少していることです。Level 1 は512点、Level 2 は256点、…、Level 4 は64点です。これはダウンサンプリングの結果であり、高周波成分ほど時間分解能が高く（点数が多く）、低周波成分ほど周波数分解能が高い（時間方向は粗い）というウェーブレット変換の適応的分解能を反映しています。

ウェーブレット変換の応用例

信号のノイズ除去（デノイジング）

ウェーブレット変換の重要な応用の一つが信号のノイズ除去です。基本的なアイデアは以下のとおりです。

信号をDWTで分解する
各レベルの詳細係数に閾値処理（thresholding）を適用する
処理後の係数から信号を再構成する

ノイズは一般に全スケールに均等に分布する一方、信号の重要な特徴は特定のスケールに集中します。したがって、小さな詳細係数（ノイズに対応）をゼロにするハード閾値処理、または縮小するソフト閾値処理により、信号成分を保持しつつノイズを除去できます。

ソフト閾値処理は

$$ \eta_\lambda(d) = \text{sgn}(d) \cdot \max(|d| – \lambda, 0) $$

で定義されます。閾値 $\lambda$ の選択には、Donoho-Johnstonのユニバーサル閾値 $\lambda = \sigma\sqrt{2\log N}$ がよく使われます。ここで $\sigma$ はノイズの標準偏差の推定値、$N$ はデータ点数です。$\sigma$ は最も高い周波数帯域（Level 1）の詳細係数のメディアン絶対偏差（MAD）から $\sigma \approx \text{MAD}/0.6745$ として推定できます。

画像圧縮（JPEG 2000）

JPEG 2000規格は2次元DWTに基づく画像圧縮方式です。従来のJPEG（DCTベース）と比較して以下の利点があります。

ブロックノイズの排除: JPEGは $8 \times 8$ ブロック単位で処理するためブロック境界にアーティファクトが出ますが、JPEG 2000は画像全体にウェーブレット変換を適用するため境界アーティファクトが生じません
プログレッシブ表示: 低解像度から高解像度へ段階的にデータを送信できます。これは多重解像度解析の自然な帰結です
高圧縮率での品質: 低ビットレートにおいて、JPEG 2000はJPEGよりも主観的な画質が良好です

2次元DWTは、行方向のDWTと列方向のDWTを組み合わせて実現します。1レベルの分解で画像は4つのサブバンド（LL, LH, HL, HH）に分かれます。LL は低周波成分（縮小画像）、LH は水平エッジ、HL は垂直エッジ、HH は対角エッジの情報を含みます。

金融データの多重スケール解析

株価や為替レートの時系列データは、短期変動（日中のノイズ）、中期トレンド（週〜月単位の動き）、長期トレンド（年単位の構造変化）が重なった複雑な信号です。ウェーブレット変換により、これらの異なるスケールの成分を分離し、それぞれを独立に分析することが可能になります。

たとえば、ウェーブレットの詳細係数のエネルギー分布が時間とともに変化するパターンを監視することで、市場のボラティリティ変動をスケールごとに捕捉できます。

まとめ

本記事では、ウェーブレット変換の基礎理論から実装、応用まで体系的に解説しました。

フーリエ変換の限界: フーリエ変換は時間情報を失い、STFTは固定窓による時間-周波数分解能のトレードオフがある。不確定性原理 $\Delta t \cdot \Delta\omega \geq 1/2$ が根本的な制約である
ウェーブレット変換はスケーリングとシフトされたマザーウェーブレットとの内積で、時間-周波数の両方の情報を同時に捉える。低周波で良い周波数分解能、高周波で良い時間分解能という適応的な分解能を持つ
CWTは連続的なスケールとシフトで完全な時間-周波数表現を与える。許容条件 $C_\psi < \infty$ を満たすゼロ平均のウェーブレットが必要
マザーウェーブレットは目的に応じて選択する。周波数分析にはMorlet、エッジ検出にはMexican hat、高速計算にはHaar/Daubechiesが適する
DWTはフィルタバンク（Mallatのアルゴリズム）として $O(N)$ で効率的に実装でき、多重解像度解析の基盤となる。完全再構成条件を満たすフィルタ設計が重要である
応用: ノイズ除去（閾値処理）、画像圧縮（JPEG 2000）、金融データの多重スケール解析、地震波解析、医療信号処理（ECG異常検出）など幅広い分野で活用されている

次のステップとして、以下の記事も参考にしてください。

フーリエ変換の定義と性質 — 周波数解析の基本
1階常微分方程式の解法まとめ — 微分方程式への展開
畳み込み定理とフーリエ変換の関係 — 畳み込みの理論