はじめに
ウェーブレット変換入門 で扱った離散ウェーブレット変換(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\) まで展開したとき DWT の葉ノード数は \(J + 1\) (各レベルの \(d\) 分岐 1 個ずつ + 最終レベルの \(a\) 分岐)であるのに対し、WPT の葉ノード数は \(2^J\) 個です。\(J = 5\) なら DWT は 6 枚、WPT は 32 枚の葉候補を持つことになり、この差が後述する「基底辞書」の大きさに直結します。
レベル \(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)}")
実行結果は次の通りです。
レベル5のノード数: 32
各ノードの係数長: 38
db4(フィルタ長 8)で symmetric 境界拡張を使うと、ダウンサンプリングのたびに端点処理分の係数が余分に生じるため、単純な \(N / 2^5 = 1024 / 32 = 32\)
よりわずかに多い 38 サンプルになります。この境界効果はレベルが深くなるほど累積するため、短い信号を深く分解する場合はこのオーバーヘッド比率が無視できなくなる点に注意が必要です。
WaveletPacket は遅延評価で、wp['aad'] のように二進パス(a=approximation, d=detail)でノードに直接アクセスできます。get_level(j, order="natural") は自然順(フィルタリング順)で \(2^j\)
個のノードを返します。
エントロピーコストの計算
シャノンエントロピーを加法的コストとして実装します。加法性を成立させる鍵は、正規化に使うエネルギーを「そのノード自身のエネルギー」ではなく「原信号全体のエネルギー \(\|x\|^2\) 」に固定することです。
def shannon_entropy_cost(coeffs, total_energy, eps=1e-12):
"""係数ベクトルのシャノンエントロピー型コスト。
p_k = |c_k|^2 / \|x\|^2 として H = -sum(p_k * log(p_k))。
total_energy には原信号全体のエネルギー sum(|x|^2) を渡す(ノード自身のエネルギーではない)。
"""
p = (np.abs(coeffs) ** 2) / total_energy
p = p[p > eps]
return float(-np.sum(p * np.log(p)))
total_energy = np.sum(np.abs(signal) ** 2)
落とし穴:ノード自身のエネルギーで正規化してはいけない
実装でもっとも起きやすい間違いは、p_k の分母に total_energy(原信号全体のエネルギー)ではなく np.sum(np.abs(coeffs) ** 2)(そのノード自身のエネルギー)を使ってしまうことです。式 (9) の直前の説明にある通り理論上は
であるべきですが、コードの見た目上はどちらも「正規化して確率にする」という自然な操作に見えるため、レビューでも見逃されがちです。
この2つがなぜ違う結果を生むかを確認します。ノード自身のエネルギーで正規化すると、各ノードの \(p_k\) は常に合計 1 になるため、
\[ \sum_k p_k \log(1/p_k) \]はそのノードの係数分布の形だけを反映し、「そのノードにどれだけのエネルギーが割り当てられているか」という情報が失われます。その結果、子ノード2つのコスト和が親ノードのコストを常に上回りやすくなり(分割してもコストが下がったように見えない)、木を深く分解する動機が消えてしまいます。実際に両方の正規化で実行すると次のようになります。
# ノード自身のエネルギーで正規化(誤り)
選択された葉ノード数: 1
最適コスト(エントロピー): 6.4792
採用されたパス例: ['']
# 原信号全体のエネルギーで正規化(正しい)
選択された葉ノード数: 17
最適コスト(エントロピー): 4.7176
採用されたパス例: ['aaaaa', 'aaaad', 'aaad', 'aad', 'adaa', 'adada', 'adadd', 'adda']
誤った正規化では Best Basis がルートノード(分解なし)を選んでしまい、ウェーブレットパケット分解を行う意味がなくなります。正しく原信号全体のエネルギーで正規化すると、過渡パルスの時間帯・低周波の正弦波・高周波の正弦波それぞれに適応した 17 個の葉ノードが選ばれます。以降のコードはすべて正しい(グローバル正規化の)実装を前提とします。
Best Basis アルゴリズムの実装
PyWavelets には Best Basis を直接返す API はないため、(10) 式の再帰を自前で実装します。
def best_basis(wp, total_energy, 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 = shannon_entropy_cost(node.data, total_energy)
# 葉(最大分解レベル)に到達した場合は自身を返す
if len(current_path) == level:
return [current_path], own_cost
# 子ノードを再帰的に評価
left_paths, left_cost = best_basis(wp, total_energy, level, current_path + "a")
right_paths, right_cost = best_basis(wp, total_energy, 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, total_energy, level=5)
print(f"選択された葉ノード数: {len(selected)}")
print(f"最適コスト(エントロピー): {total_cost:.4f}")
print("採用されたパス例:", selected[:8])
実行結果は次の通りです(上記の「正しい」ケースと同じ)。
選択された葉ノード数: 17
最適コスト(エントロピー): 4.7176
採用されたパス例: ['aaaaa', 'aaaad', 'aaad', 'aad', 'adaa', 'adada', 'adadd', 'adda']
完全分解(\(2^5 = 32\)
葉)と比較して、選択された葉ノード数は 17 と半分程度に減っています。選ばれたパスの周波数帯域を計算すると、da(レベル2、384–512Hz 帯)だけが唯一の浅い葉(時間分解能を優先)で、パルス由来の広帯域ノイズが支配的な高域を粗い周波数分解能でまとめて捉えています。一方、20Hz 正弦波を含む aaaad(レベル5、16–32Hz 帯)や 200Hz 正弦波を含む adaa(レベル4、192–224Hz 帯)は深い葉として選ばれ、狭帯域の定常成分を高い周波数分解能で分離しています。信号の構造(過渡パルス vs 定常正弦波)に応じたハイブリッドな分解になっていることが確認できます。
Best Basis の時間周波数タイル可視化
選ばれた葉の集合を時間周波数平面上のタイルとして描画します。ここで一つ注意が必要です。各葉ノードの係数列 wp[path].data は、原信号 \(N=1024\)
サンプルをレベル \(j\)
回ダウンサンプリングした \(M \approx N / 2^j\)
個の時系列であり、係数 1 個 1 個が時間軸上の異なる位置に対応します。したがって「葉 1 個 = 時間軸全体を覆う 1 枚の矩形」という単純化した描画では、周波数分解能(矩形の高さ)しか可視化できず、WPT の核心である時間分解能の違いが図から消えてしまいます。正しくは、各葉の \(M\)
個の係数それぞれを時間軸上の小タイルとして並べる必要があります。
def plot_best_basis_tf_grid(selected_paths, wp, fs, ax):
"""Best Basisの葉を、係数1個=1タイルとして時間周波数平面に描画。
タイルの濃淡は係数の絶対値(そのタイルにどれだけエネルギーが
集中しているか)を表す。
"""
for path in selected_paths:
j = len(path) # 分解レベル = 深さ
node = wp[path] if path else wp
coeffs = np.abs(node.data)
M = len(coeffs) # このノードの係数数(時間分割数)
# パスを自然順のサブバンド番号 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)
cmax = coeffs.max() if coeffs.max() > 0 else 1.0
for k in range(M):
alpha = 0.12 + 0.75 * (coeffs[k] / cmax)
ax.add_patch(
plt.Rectangle(
(k / M, f_low),
1.0 / M,
f_high - f_low,
facecolor="C0",
alpha=min(alpha, 0.87),
edgecolor="white",
linewidth=0.15,
)
)
# 帯域の輪郭線
ax.add_patch(
plt.Rectangle((0, f_low), 1.0, f_high - f_low, fill=False, edgecolor="C0", linewidth=0.8)
)
ax.set_xlim(0, 1)
ax.set_ylim(0, fs / 2)
ax.set_xlabel("Time (normalized)")
ax.set_ylabel("Frequency [Hz]")
ax.set_title(f"Best Basis time-frequency tiles ({len(selected_paths)} leaves)")
fig, ax = plt.subplots(figsize=(9, 5.8))
plot_best_basis_tf_grid(selected, wp, fs, ax)
plt.tight_layout()
plt.show()
![Best Basisの時間周波数タイル可視化。横軸は正規化時間、縦軸は周波数[Hz]、17個の葉ノードそれぞれが係数1個=1タイルとして描画され、濃淡はその係数の絶対値を表す。384〜512Hz帯(daノード、レベル2で最も浅い=時間分解能優先)にt≈0.29付近で濃い縦縞が現れており、これはパルス(サンプル300〜305、t≈0.293〜0.298)の位置と一致する。また120〜250Hz帯の複数バンドではt=0.5を境に濃淡が急変しており、t>0.5でオンになる200Hz正弦波の開始点を捉えている](/posts/20260522_wavelet_packet/best_basis_tiling.png)
図から、da(384–512Hz、深さ2)の帯では \(t \approx 0.29\)
付近に濃い縦縞が現れており、これはパルスの実際の位置(サンプル 300–305、\(t = 300/1024 \approx 0.293\)
から \(305/1024 \approx 0.298\)
)と一致します。浅いノードは時間分割数 \(M\)
が少ない(時間分解能は粗いが検出は速い)にもかかわらず、パルスのような広帯域・短時間の成分を効率よく捉えられることが分かります。また 120–250Hz 帯の複数バンドでは \(t = 0.5\)
を境に濃淡が明確に変化しており、t > 0.5 でオンになる 200Hz 正弦波の開始点を正しく反映しています。これは
STFT
の固定タイル(全タイルが同じ形)や DWT の左偏ったタイル分割にはない、帯域ごとに時間分解能と周波数分解能の配分を変えるという WPT Best Basis の本質を示しています。
完全分解 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が最小化を保証する量)
level5_cost = sum(shannon_entropy_cost(n.data, total_energy) for n in nodes_level5)
print(f"レベル5固定分解の合計エントロピー: {level5_cost:.4f}")
print(f"Best Basisの合計エントロピー: {total_cost:.4f}")
実行結果は次の通りです。
完全分解(レベル5固定): 上位5%でエネルギー 0.846
Best Basis : 上位5%でエネルギー 0.848
レベル5固定分解の合計エントロピー: 4.8372
Best Basisの合計エントロピー: 4.7176
ここで注意すべき点があります。上位 5% エネルギー集中度は 0.846 → 0.848 とわずかしか改善していません。実際、\(k\)
の値を 1%・2%・10% に変えて計測すると、Best Basis の方がむしろ低い場合すらあります(1%: 完全分解 0.396 / Best Basis 0.380、2%: 0.613 / 0.592、10%: 0.949 / 0.941)。これは実装のバグではなく、「上位 \(k\)
% エネルギー集中度」と「加法的エントロピーコスト」は異なる量を測る指標であることに起因します。Best Basis アルゴリズムが数学的に保証するのは、加法的コスト関数 \(\mathcal{M}(\mathcal{T})\)
(ここではエントロピー)の最小化であって、上位 \(k\)
% エネルギー集中度のような別の指標の最大化ではありません。実際、合計エントロピーで比較すると、完全分解(レベル5固定)の 4.8372 に対して Best Basis は 4.7176 と、Best Basis の構成上の保証通り確実に下回っています。加えて、mags_full と mags_best は要素数が異なる(それぞれ 1216 個・1125 個)ため、上位 \(k\)
% が指す絶対個数もずれており、この種の相対指標を異なる基底同士で比較する際は絶対個数を揃えるなど注意が必要です。スパース性を評価する際は、Best Basis が実際に最小化している目的関数(エントロピーコスト)そのもので比較するのが最も信頼できる方法です。
用途
ウェーブレットパケット変換と Best Basis 選択は、次のような場面で特に有効です。
信号圧縮
JPEG 2000 などで使われる DWT ベースの圧縮を一般化し、信号ごとに最適な基底でスパース表現する。狭帯域の周波数構造が支配的な音響信号などで、DWT より高い圧縮率を達成できる場合があります。
非定常信号の分類
EEG(脳波)、EMG(筋電図)、機械振動などでは、各クラスに特徴的な周波数帯域・時間帯域が存在します。Best Basis でクラスごとに代表的なパケットを選び、その係数統計を特徴量とすると、固定タイル分割では捉えにくいクラス間差異を抽出できます。これは 時系列データの異常検知 の特徴抽出としても応用されます。
音声・音響解析
音声は声道共鳴に対応する狭帯域構造(フォルマント)と過渡的な破裂音を併せ持つため、固定窓の STFT では一方を犠牲にせざるを得ません。Best Basis は両方を同時に高分解能で捉えることができ、音楽の和音解析やイルカ・コウモリの非定常音の分析にも使われます。
通信・レーダー信号の検出
未知の変調方式・チャープ波形の検出では、エネルギー集中度の高い葉ノードを特徴として使うことで、 FFT や 窓関数・PSD だけでは見えない構造を抽出できます。
比較表:DWT vs WPT + Best Basis
| 特性 | DWT | WPT + 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 に対して常に等しいかそれ以上の表現効率を持つ保証があります。
最新研究動向:非デシメーション型ウェーブレットパケット特徴量と深層学習予測
ウェーブレットパケット変換は古典的な信号圧縮・分類にとどまらず、近年は深層学習ベースの時系列予測モデルの入力特徴量としても再評価されています。
Nason & Wei (2024)「 Leveraging Non-Decimated Wavelet Packet Features and Transformer Models for Time Series Forecasting 」(arXiv:2403.08630)は、通常の(間引きを伴う)ウェーブレットパケット変換ではなく、ダウンサンプリングを行わない非デシメーション型ウェーブレットパケット変換の係数を特徴量として用い、統計的手法から Transformer ベースの深層学習モデルまで幅広い予測手法で評価しています。情報漏洩(未来の情報が特徴量に混入すること)を防ぐため、シフトされたピラミッド型アルゴリズムを用いて各時刻の特徴量を計算する点が工夫です。実験では、非時系列的な(各時刻を独立に扱う)予測手法において、従来の高次ラグ特徴量をウェーブレットパケット特徴量に置き換えることで単一ステップ予測に大きな改善が見られた一方、Transformer などの時系列深層学習モデルではより長い予測期間に対する改善は穏やかだったと報告されています。本記事で解説した「完全二分木による全帯域分解」という WPT の構造が、間引きなしの形でも特徴量エンジニアリングに応用できる点は、Best Basis 選択とは異なるもう一つの実用的な活用方向として興味深い発展です。
まとめ
- ウェーブレットパケット変換は DWT が分解しない高周波側も再帰的に二分割し、\(2^J\) 個のサブ空間に直交分解する。
- 完全分解を成す葉ノード集合を 1 つ選ぶごとに正規直交基底が定まり、ツリー全体は基底辞書とみなせる。
- Coifman–Wickerhauser の Best Basis アルゴリズムは加法的コスト関数(シャノンエントロピー等)を \(O(N \log N)\) で最小化し、信号適応的な最適基底を選び出す。ただし正規化は原信号全体のエネルギーで行う必要があり、ノード自身のエネルギーで正規化すると分解が機能しなくなる(実測: 誤りでは葉1個・コスト6.4792、正しくは葉17個・コスト4.7176)。
- PyWavelets の
WaveletPacketと自前の再帰実装を組み合わせることで、Best Basis 選択・タイル可視化・スパース性評価を一貫して実装できる。ただし Best Basis が保証するのはエントロピーコストの最小化であり、上位 \(k\) % エネルギー集中度のような別の指標は必ずしも改善するとは限らない。 - 用途は信号圧縮、非定常信号の分類、音声・通信信号解析と幅広く、DWT のサブセットとして常に等しいかそれ以上の効率を持つ。
ウェーブレットパケット変換と Best Basis は、ウェーブレット理論を信号適応的な時間周波数解析へと拡張する自然な一歩であり、固定タイル分割の限界を超える鍵となる枠組みです。
関連記事
- ウェーブレット変換入門:時間-周波数解析のPython実装 - 本記事の前提となる CWT/DWT・MRA・フィルタバンクを解説しています。
- 短時間フーリエ変換(STFT)の原理とPython実装 - 固定窓による時間周波数解析。WPT + Best Basis が信号適応的にタイル形を変える比較対象となります。
- 高速フーリエ変換(FFT)の仕組みとPython実装 - 周波数解析の基礎であり、各パケットノード内のスペクトル特性を考える出発点です。
- DTFT・DFT・FFTの違いを整理する - 離散信号の周波数表現の階層を整理しており、パケット係数の周波数解釈と接続します。
- 窓関数とパワースペクトル密度(PSD)の理論とPython実装 - パケットノード内でのスペクトル推定の前提となる窓関数と PSD を解説しています。
- ヒルベルト変換と解析信号:瞬時振幅・位相・周波数のPython実装 - 各パケット係数の包絡・瞬時周波数を取り出す解析ツールとして併用できます。
- 時系列データの異常検知 - Best Basis 係数を特徴量とした特徴抽出は異常検知の有力な前処理になります。
- 時間周波数解析の選び方ハブ - WPT を FFT・STFT・CWT・Hilbert変換と並べて、タイル形と信号適応性の観点で位置付ける統合ハブです。
- MDCT(修正離散コサイン変換)とフィルタバンク - WPT の2分割フィルタバンクの木構造に対し、M帯域を一括で均等分割するコサイン変調フィルタバンクの代表例です。
参考文献
- 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.
- Nason, G. P., & Wei, J. L. (2024). Leveraging Non-Decimated Wavelet Packet Features and Transformer Models for Time Series Forecasting . arXiv:2403.08630.
- PyWavelets documentation: WaveletPacket