ウェーブレットパケット変換と最適基底選択：時間周波数解析の高分解能化とPython実装

ウェーブレットパケット変換の基礎・設計・Python実装とDWTとの比較を解説。Best Basis選択による高分解能な時間周波数解析を、Coifman-Wickerhauserエントロピーとバイナリツリー分解でPyWaveletsで実装します。

5 min read

はじめに

ウェーブレット変換入門 で扱った離散ウェーブレット変換（DWT）は、信号を「近似（低周波）」と「詳細（高周波）」に再帰的に分解する強力な手法です。しかし、DWT の標準アルゴリズム（Mallat 分解）には次のような制約があります。

DWT は近似側（低周波側）だけを再帰的に分解し、詳細側（高周波側）は一度しか分解しない。このため、高周波帯域の周波数分解能が粗く、たとえば 8kHz サンプリングの信号を 4 段分解した場合、最終レベルの詳細係数（cD1）は 2–4 kHz 帯域に対応しますが、この広い帯域の中での周波数構造は捉えられません。

ウェーブレットパケット変換（Wavelet Packet Transform, WPT）はこの制約を取り払い、詳細側も含めて全帯域を再帰的に二分割することで、時間周波数平面を柔軟かつ高分解能に分解します。さらに、コスト関数（典型的にはエントロピー）に基づいて最良のサブツリー＝**最適基底（Best Basis）**を選び出すことで、信号に応じた最適な時間周波数分解を自動的に決定できます。

本記事では、ウェーブレットパケットの数学的構造、Coifman–Wickerhauser の Best Basis アルゴリズム、そして PyWavelets を用いた Python 実装を解説します。

ウェーブレットパケットの基本

バイナリツリーによる完全分解

DWT が「左偏った木」（低周波側のみを再帰分解）であるのに対し、ウェーブレットパケット変換は完全二分木を成します。各ノードはローパスフィルタ \(h\) とハイパスフィルタ \(g\) により 2 つの子ノードに分解されます。

レベル \(j\) のノード \((j, n)\) のパケット係数を \(d_{j,n}[k]\) とすると、次の漸化式が成り立ちます。

\[d_{j+1, 2n}[k] = \sum_{l} h[l - 2k] \, d_{j, n}[l] \tag{1}\] \[d_{j+1, 2n+1}[k] = \sum_{l} g[l - 2k] \, d_{j, n}[l] \tag{2}\]

ここで \(n\) は同一レベル内のサブバンド番号（\(0 \le n < 2^j\) ）、\(h\) ・\(g\) はそれぞれスケーリングフィルタ・ウェーブレットフィルタです。ルート \((0, 0)\) は原信号 \(x[k]\) に対応します。

ウェーブレットパケット基底関数

連続時間で見ると、ウェーブレットパケット基底 \(\psi_{n,j,k}(t)\) は次の二重スケール方程式から再帰的に生成されます。

\[\psi_{2n}(t) = \sqrt{2} \sum_k h[k] \, \psi_n(2t - k) \tag{3}\] \[\psi_{2n+1}(t) = \sqrt{2} \sum_k g[k] \, \psi_n(2t - k) \tag{4}\]

ここで \(\psi_0 = \phi\) （スケーリング関数）、\(\psi_1 = \psi\) （マザーウェーブレット）です。シフトとダイアディックスケーリングを加えた基底 \(\psi_{n,j,k}(t) = 2^{-j/2}\psi_n(2^{-j} t - k)\) は、レベル \(j\) のサブバンド \(n\) における基底関数を表します。

DWT との関係は単純で、DWT は各レベルで \(n = 0\) （近似側）だけを分解する特殊な部分木に相当します。WPT はこの制限を外し、全ノードを分解する完全二分木を扱います。

完全分解と冗長性

サブ空間の階層構造

完全に分解レベル \(J\) まで展開すると、\(L^2(\mathbb{R})\) 空間は \(2^J\) 個の直交サブ空間 \(W_{J,n}\) （\(0 \le n < 2^J\) ）に分解されます。

\[L^2(\mathbb{R}) = \bigoplus_{n=0}^{2^J - 1} W_{J,n} \tag{5}\]

各 \(W_{J,n}\) は周波数軸を等分割した帯域に対応し、レベル \(J\) のフィルタバンク出力 \(d_{J,n}\) がその係数となります。同時に、サブ空間には 直交分解の入れ子が成り立ちます。

\[W_{j,n} = W_{j+1, 2n} \oplus W_{j+1, 2n+1} \tag{6}\]

つまり、ある親ノードの基底を「子ノード 2 つの基底に置き換える」ことが直交変換として正当化されます。これがあとで Best Basis アルゴリズムの根拠になります。

エネルギー保存則（Parseval）

ウェーブレットパケット変換は直交変換であるため、各サブツリーが完全分解（葉ノード集合が \([0, N)\) を被覆）である限り、Parseval の等式によって信号エネルギーが保存されます。

\[\sum_{k} |x[k]|^2 = \sum_{(j, n) \in \mathcal{T}} \sum_{k} |d_{j,n}[k]|^2 \tag{7}\]

ここで \(\mathcal{T}\) は完全分解を成す葉ノードの集合です。エネルギー保存則は、エントロピー型コストの比較を意味あるものにする鍵となる性質です（後述）。

冗長性と「基底」の選び方

ツリー全体には \((2J + 1) \cdot N\) 個程度の係数があり、原信号の \(N\) サンプルに対して冗長です。しかし、完全分解を成す葉ノード集合（前述の \(\mathcal{T}\) ）を 1 つ選べば、その集合の係数だけで原信号を表現できる正規直交基底になる点が WPT の核心です。

つまり WPT のツリーは「無数の正規直交基底の辞書」を提供し、その中から信号に最も適した 1 つを選ぶのが次節の Best Basis 選択です。

最適基底選択（Best Basis）

コスト関数とエントロピー

信号 \(x\) をある基底（葉ノード集合 \(\mathcal{T}\) ）で表現したときの「複雑さ」を測るスカラー量を加法的コスト関数 \(\mathcal{M}(\mathcal{T})\) と呼びます。加法的とは、\(\mathcal{T}\) に属するすべての係数 \(d_{j,n}[k]\) について

\[\mathcal{M}(\mathcal{T}) = \sum_{(j, n) \in \mathcal{T}} \sum_{k} \mu(|d_{j,n}[k]|^2) \tag{8}\]

の形で書けることを意味します。代表的なコスト関数として、Coifman と Wickerhauser が導入したシャノンエントロピー型コストがあります。

\[\mu(p) = -p \log p \quad (p > 0), \quad \mu(0) = 0 \tag{9}\]

ここで \(p_k = |d_{j,n}[k]|^2 / \|x\|^2\) と正規化すれば、\(\mathcal{M}\) は係数分布の不確定性を表すシャノン情報エントロピーとなり、値が小さいほど「少数の大きな係数にエネルギーが集中している＝スパース」を意味します。

他にも、\(\ell^p\) ノルム（\(p < 2\) ）、対数エネルギー \(\sum \log |d|^2\) 、閾値超え個数などが使われます。共通の要件は加法性で、これにより部分木ごとのコスト比較が局所的に完結します。

Best Basis アルゴリズム

Coifman–Wickerhauser の Best Basis アルゴリズム（1992）は、すべての完全分解 \(\mathcal{T}\) の中で \(\mathcal{M}(\mathcal{T})\) を最小化するものを \(O(N \log N)\) で見つけます。ボトムアップの動的計画法で、各ノードについて次の決定を再帰的に行います。

\[\mathcal{M}_{\text{best}}(j, n) = \min \big\{ \mathcal{M}(d_{j, n}), \, \mathcal{M}_{\text{best}}(j+1, 2n) + \mathcal{M}_{\text{best}}(j+1, 2n+1) \big\} \tag{10}\]
  • 左辺：そのノードを葉として採用したときのコスト
  • 右辺第 2 項：子 2 つの最適コストの和

子の合計コストが親自身のコストを下回るならば「子に分解する」、そうでなければ「このノードを葉として残す」と決めます。葉ノードでは比較せず、自身のコストを返します。葉から根へ向かって 1 回走査すれば最適基底が決まるため、計算量は係数生成と同じ \(O(N \log N)\) です。

Best Basis の幾何的意味

エントロピー最小化は、信号エネルギーが「少数のパケット係数に集中する」分解を選ぶことに相当します。

  • 純粋な正弦波 → 周波数局在が良い狭帯域パケットに集中 → 周波数解像度の高い葉が選ばれる
  • インパルス・過渡波形 → 時間局在が良い浅い葉が選ばれる
  • チャープ・周波数変化信号 → 時間と周波数の中間的な葉が選ばれる

つまり Best Basis は、STFT のような固定タイル分割や DWT のような左偏ったタイル分割を超え、信号適応的に時間周波数タイルの形を選ぶ仕組みです。

Python 実装

PyWavelets による完全分解

pywt.WaveletPacket を使うと、完全分解と任意ノードへのアクセスが簡単に行えます。

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

# --- テスト信号: 2つの正弦波 + 過渡パルス ---
np.random.seed(0)
fs = 1024
N = 1024
t = np.arange(N) / fs

low = np.sin(2 * np.pi * 20 * t)
high = 0.6 * np.sin(2 * np.pi * 200 * t) * (t > 0.5)
pulse = np.zeros(N)
pulse[300:305] = 2.0
signal = low + high + pulse + 0.05 * np.random.randn(N)

# --- ウェーブレットパケットの完全分解 ---
wp = pywt.WaveletPacket(data=signal, wavelet="db4", mode="symmetric", maxlevel=5)

# レベル5の全ノードを取得（Naturalオーダー）
nodes_level5 = wp.get_level(5, order="natural")
print(f"レベル5のノード数: {len(nodes_level5)}")  # 2^5 = 32
print(f"各ノードの係数長: {len(nodes_level5[0].data)}")

WaveletPacket は遅延評価で、wp['aad'] のように二進パス（a=approximation, d=detail）でノードに直接アクセスできます。get_level(j, order="natural") は自然順（フィルタリング順）で \(2^j\) 個のノードを返します。

エントロピーコストの計算

シャノンエントロピーを加法的コストとして実装します。

def shannon_entropy_cost(coeffs, eps=1e-12):
    """係数ベクトルのシャノンエントロピー型コスト。

    p_k = |c_k|^2 / sum(|c|^2) として H = -sum(p_k * log(p_k))。
    エネルギー保存則のもとで加法的コストになる。
    """
    energy = np.sum(np.abs(coeffs) ** 2)
    if energy < eps:
        return 0.0
    p = (np.abs(coeffs) ** 2) / energy
    p = p[p > eps]
    return float(-np.sum(p * np.log(p)))

Best Basis アルゴリズムの実装

PyWavelets には Best Basis を直接返す API はないため、(10) 式の再帰を自前で実装します。

def best_basis(wp, cost_fn, level, current_path=""):
    """Coifman-Wickerhauser Best Basis を再帰的に探索する。

    Returns
    -------
    selected_paths : list[str]   採用された葉ノードのパス
    total_cost     : float       そのサブツリーの最適コスト
    """
    node = wp[current_path] if current_path else wp
    own_cost = cost_fn(node.data)

    # 葉（最大分解レベル）に到達した場合は自身を返す
    if len(current_path) == level:
        return [current_path], own_cost

    # 子ノードを再帰的に評価
    left_paths, left_cost = best_basis(wp, cost_fn, level, current_path + "a")
    right_paths, right_cost = best_basis(wp, cost_fn, level, current_path + "d")
    child_cost = left_cost + right_cost

    if child_cost < own_cost:
        return left_paths + right_paths, child_cost
    else:
        return [current_path], own_cost


selected, total_cost = best_basis(wp, shannon_entropy_cost, level=5)
print(f"選択された葉ノード数: {len(selected)}")
print(f"最適コスト（エントロピー）: {total_cost:.4f}")
print("採用されたパス例:", selected[:8])

完全分解（\(2^5 = 32\) 葉）と比較して、選択された葉ノード数が少なくなるはずです。信号の構造によっては、過渡パルスが存在する時間帯では浅いノード、定常正弦波が支配的な帯域では深いノードが選ばれる、というハイブリッドな分解が得られます。

Best Basis の時間周波数タイル可視化

選ばれた葉の集合を時間周波数平面上のタイルとして描画します。

def plot_best_basis_tiling(selected_paths, max_level, fs, ax=None):
    """Best Basisの葉を時間周波数タイルとして可視化。"""
    if ax is None:
        _, ax = plt.subplots(figsize=(10, 5))

    for path in selected_paths:
        j = len(path)  # 分解レベル = 深さ
        # パスをグレイコード順のサブバンド番号nに変換（自然順を周波数順に並べる）
        n_natural = int(path.replace("a", "0").replace("d", "1"), 2) if path else 0
        # 自然順 -> 周波数順（グレイコード逆変換）
        n_freq = n_natural
        for shift in (1,):
            n_freq ^= n_natural >> shift

        f_low = n_freq * (fs / 2) / (2 ** j)
        f_high = (n_freq + 1) * (fs / 2) / (2 ** j)
        ax.add_patch(
            plt.Rectangle(
                (0, f_low),
                1.0,  # 葉は時間全域を覆う（パケット係数列を時間軸として展開する場合は別途）
                f_high - f_low,
                fill=False,
                edgecolor="C0",
                linewidth=1.0,
            )
        )

    ax.set_xlim(0, 1)
    ax.set_ylim(0, fs / 2)
    ax.set_xlabel("Time (normalized)")
    ax.set_ylabel("Frequency [Hz]")
    ax.set_title("Best Basis Tiling (Wavelet Packet)")
    return ax


fig, ax = plt.subplots(figsize=(10, 5))
plot_best_basis_tiling(selected, max_level=5, fs=fs, ax=ax)
plt.tight_layout()
plt.show()

タイルは「縦長 = 周波数分解能を犠牲にして時間分解能を取った」「横長 = 逆」を意味し、信号構造に応じて適応的に変化します。これは STFT の固定タイル（全タイルが同じアスペクト比）や DWT の左偏ったタイル分割と本質的に異なる点です。

完全分解 vs Best Basis のエネルギー集中度

最終的に Best Basis がスパース表現として優れているかを定量評価します。

def coefficient_magnitudes(wp, paths):
    """選択されたパス集合の全係数を1次元配列にまとめる。"""
    coeffs = []
    for p in paths:
        coeffs.append(np.abs(wp[p].data))
    return np.concatenate(coeffs)


full_paths = [n.path for n in wp.get_level(5, order="natural")]
mags_full = np.sort(coefficient_magnitudes(wp, full_paths))[::-1]
mags_best = np.sort(coefficient_magnitudes(wp, selected))[::-1]


def energy_concentration(mags, k_ratio=0.05):
    """上位 k_ratio の係数が捉えるエネルギー比。"""
    k = max(1, int(len(mags) * k_ratio))
    total = np.sum(mags ** 2)
    top = np.sum(mags[:k] ** 2)
    return top / total


print(f"完全分解（レベル5固定）: 上位5%でエネルギー {energy_concentration(mags_full):.3f}")
print(f"Best Basis            : 上位5%でエネルギー {energy_concentration(mags_best):.3f}")

Best Basis は信号に適応してエネルギーを少数の係数に集中させるため、圧縮や閾値処理におけるスパース性の指標が改善します。

用途

ウェーブレットパケット変換と Best Basis 選択は、次のような場面で特に有効です。

信号圧縮

JPEG 2000 などで使われる DWT ベースの圧縮を一般化し、信号ごとに最適な基底でスパース表現する。狭帯域の周波数構造が支配的な音響信号などで、DWT より高い圧縮率を達成できる場合があります。

非定常信号の分類

EEG（脳波）、EMG（筋電図）、機械振動などでは、各クラスに特徴的な周波数帯域・時間帯域が存在します。Best Basis でクラスごとに代表的なパケットを選び、その係数統計を特徴量とすると、固定タイル分割では捉えにくいクラス間差異を抽出できます。これは 時系列データの異常検知 の特徴抽出としても応用されます。

音声・音響解析

音声は声道共鳴に対応する狭帯域構造（フォルマント）と過渡的な破裂音を併せ持つため、固定窓の STFT では一方を犠牲にせざるを得ません。Best Basis は両方を同時に高分解能で捉えることができ、音楽の和音解析やイルカ・コウモリの非定常音の分析にも使われます。

通信・レーダー信号の検出

未知の変調方式・チャープ波形の検出では、エネルギー集中度の高い葉ノードを特徴として使うことで、FFT や 窓関数・PSD だけでは見えない構造を抽出できます。

比較表：DWT vs WPT + Best Basis

特性DWTWPT + Best Basis
分解構造左偏った木（近似側のみ分解）完全二分木（全帯域を再帰分解）
高周波分解能粗い（最終 cD1 が広帯域）細かい（任意レベルまで分割可能）
基底固定信号適応的（コスト最小化で自動選択）
基底候補数1\(\sim 2^{2^{J-1}}\) 個の完全分解
計算量（分解）\(O(N)\)\(O(N \log N)\)
計算量（基底）不要\(O(N \log N)\) （Coifman–Wickerhauser）
エネルギー保存ありあり（葉ノード集合が完全分解の場合）
主な用途ノイズ除去・MRA圧縮・分類・スパース表現

WPT は DWT の真の拡張であり、\(J\) を十分大きく取れば DWT の分解は必ず WPT の選択肢の中に含まれます。Best Basis アルゴリズムは「DWT が最適である信号」では DWT 型の左偏った木を選ぶため、WPT は DWT に対して常に等しいかそれ以上の表現効率を持つ保証があります。

まとめ

  • ウェーブレットパケット変換は DWT が分解しない高周波側も再帰的に二分割し、\(2^J\) 個のサブ空間に直交分解する。
  • 完全分解を成す葉ノード集合を 1 つ選ぶごとに正規直交基底が定まり、ツリー全体は基底辞書とみなせる。
  • Coifman–Wickerhauser の Best Basis アルゴリズムは加法的コスト関数（シャノンエントロピー等）を \(O(N \log N)\) で最小化し、信号適応的な最適基底を選び出す。
  • PyWavelets の WaveletPacket と自前の再帰実装を組み合わせることで、Best Basis 選択・タイル可視化・スパース性評価を一貫して実装できる。
  • 用途は信号圧縮、非定常信号の分類、音声・通信信号解析と幅広く、DWT のサブセットとして常に等しいかそれ以上の効率を持つ。

ウェーブレットパケット変換と Best Basis は、ウェーブレット理論を信号適応的な時間周波数解析へと拡張する自然な一歩であり、固定タイル分割の限界を超える鍵となる枠組みです。

関連記事

参考文献

  • Coifman, R. R., & Wickerhauser, M. V. (1992). Entropy-based algorithms for best basis selection. IEEE Transactions on Information Theory, 38(2), 713–718.
  • Mallat, S. (2008). A Wavelet Tour of Signal Processing: The Sparse Way (3rd ed.). Academic Press.
  • Wickerhauser, M. V. (1994). Adapted Wavelet Analysis from Theory to Software. A K Peters.
  • PyWavelets documentation: WaveletPacket

この記事に関連する無料ツール

DevToolBox Moving Average Calculator SMA/EMA/WMA をブラウザで即座に計算 DevToolBox Frequency Converter Hz/kHz/MHz/rad/s 単位を相互変換 DevToolBox JSON Formatter JSONを整形・圧縮・バリデーション

関連記事