ウェーブレット パケット

ウェーブレット パケットの手法はウェーブレット分解の一般化であり、より豊かな信号解析を提供します。

ウェーブレット パケットの原子は、位置、(ウェーブレット分解と同様の) スケール、周波数の 3 つの自然に解釈されるパラメーターでインデックス付けされた波形です。

任意の直交ウェーブレット関数に対し、"ウェーブレット パケット基底" と呼ばれる基底のライブラリを生成します。これらの基底はそれぞれ、信号を符号化し、グローバル エネルギーを保存し、正確な特徴を再構成する特定の方法を提供します。ウェーブレット パケットを使用して、与えられた信号をさまざまに展開できます。その後、エントロピーベースの基準で、与えられた信号の最適な分解を選択します。

ウェーブレット パケット分解と最適な分解の選択の両方にシンプルで効率的なアルゴリズムが存在します。さらに、最適な信号符号化およびデータ圧縮での直接の応用で、適応フィルター処理アルゴリズムを生成できます。

ウェーブレットからウェーブレット パケットへ

直交ウェーブレット分解の処理において、一般的な手順では Approximation 係数を 2 つの部分に分割します。分割後、Approximation 係数から成るベクトルと Detail 係数から成るベクトルが、ともに粗いスケールで得られます。連続する 2 つの Approximation の間で失われた情報は、Detail 係数に取得されています。さらに、次の手順は新しい Approximation 係数ベクトルの分割で構成されます。一連の Detail が再解析されることはありません。

ウェーブレット パケットにおいても同様に、Approximation ベクトルの分割と同じ方法を使用して、各 Detail 係数ベクトルも 2 つの部分に分解されます。これにより、極めて豊かな解析が提供されます。次の図に示すような完全二分木が生成されます。

レベル 3 でのウェーブレット パケット分解ツリー

この分解の考え方は、スケール指向の分解から始めて、得られた信号を周波数サブバンドで解析することです。

実際のウェーブレット パケット: はじめに

以下の簡単な例で、ウェーブレット解析とウェーブレット パケット解析の違いを説明します。

ウェーブレット パケット スペクトル

フーリエ変換を使用する広義定常信号のスペクトル解析は、十分に確立されています。非定常信号については、短時間フーリエ変換 (STFT) などの局所フーリエ手法が存在します。概要については、短時間フーリエ変換を参照してください。

ウェーブレットは時間および周波数で局所化されるので、STFT に相当するウェーブレットベースのものを非定常信号の時間-周波数解析に使用することができます。たとえば、連続ウェーブレット変換 (CWT) に基づいたスカログラム (wscalogram) を構築できます。しかし、CWT を使用する潜在的な欠点は、計算に時間がかかることです。

離散ウェーブレット変換 (DWT) は入力信号の時間-周波数分解ができますが、DWT での周波数分解能の度合いは、一般に実用的な時間-周波数解析には粗すぎると考えられています。

DWT ベースの手法と CWT ベースの手法の間の折衷案として、ウェーブレット パケットは、計算効率が高く、十分な周波数分解能がある代替手法を提供します。wpspectrum を使用して、ウェーブレット パケットを使用した信号の時間-周波数解析を実行できます。

以下の例では、ウェーブレット パケットを使用して局所スペクトル解析を実行する方法について説明します。また、以下の例では、ウェーブレット パケット スペクトルと比較するベンチマークとして、Signal Processing Toolbox™ ソフトウェアの spectrogram (Signal Processing Toolbox) も使用します。Signal Processing Toolbox ソフトウェアがない場合は、ウェーブレット パケット スペクトルの例のみを実行します。

正弦波のウェーブレット パケット スペクトル。

fs = 1000; % sampling rate
t = 0:1/fs:2; % 2 secs at 1kHz sample rate
y = sin(256*pi*t); % sine of period 128
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

Signal Processing Toolbox ソフトウェアがある場合、短時間フーリエ変換を計算できます。

figure;
windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
[S,F,T] = spectrogram(y,window,noverlap,nfft,fs);
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

周波数が 64 Hz と 128 Hz の 2 つの正弦波の和。

fs = 1000;
t = 0:1/fs:2;
y = sin(128*pi*t) + sin(256*pi*t); % sine of periods 64 and 128.
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

Signal Processing Toolbox ソフトウェアがある場合、短時間フーリエ変換を計算できます。

figure;
windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
[S,F,T] = spectrogram(y,window,noverlap,nfft,fs);
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

16 Hz から 64 Hz への急激な周波数の変化が 2 秒の時点にある信号。

fs = 500;
t = 0:1/fs:4;
y = sin(32*pi*t).*(t<2) + sin(128*pi*t).*(t>=2);
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

Signal Processing Toolbox ソフトウェアがある場合、短時間フーリエ変換を計算できます。

figure;
windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
[S,F,T] = spectrogram(y,window,noverlap,nfft,fs);
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

線形チャープのウェーブレット パケット スペクトル。

fs = 1000;
t = 0:1/fs:2;
y = sin(256*pi*t.^2);
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

Signal Processing Toolbox ソフトウェアがある場合、短時間フーリエ変換を計算できます。

figure;
windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
[S,F,T] = spectrogram(y,window,noverlap,nfft,fs);
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

二次チャープのウェーブレット パケット スペクトル。

y = wnoise('quadchirp',10);
len = length(y);
t = linspace(0,5,len);
fs = 1/t(2);
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

Signal Processing Toolbox ソフトウェアがある場合、短時間フーリエ変換を計算できます。

windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

ウェーブレット パケットの作成

ウェーブレット パケット生成の計算手法は、直交ウェーブレットを使用すると簡単です。ウェーブレットに対応する長さ 2N の 2 つのフィルター h(n) と g(n) から始めます。

帰納法により、次の関数列を定義します。

(W_n(x), n = 0, 1, 2, ...)

定義:

$$W_{2n}(x)=\sqrt{2}{\displaystyle \sum _{k=0}^{2N-1}h(k)W_n(2x-k)}$$

$$W_{2n+1}(x)=\sqrt{2}{\displaystyle \sum _{k=0}^{2N-1}g(k)W_n(2x-k)}$$

ここで、W₀(x) = φ(x) はスケーリング関数、W₁(x) = ψ(x) はウェーブレット関数です。

たとえば、Haar ウェーブレットの場合、次が成り立ちます。

$$N=1,h(0)=h(1)=\frac{1}{\sqrt{2}}$$

と

$$g(0)=-g(1)=\frac{1}{\sqrt{2}}$$

方程式は次のようになります。

$$W_{2n}(x)=W_n(2x)+W_n(2x-1)$$

と

$$W_{2n+1}(x)=W_n(2x)-W_n(2x-1)$$

W₀(x) = φ(x) は Haar スケーリング関数で、W₁(x) = ψ(x) は Haar ウェーブレットです。どちらも [0, 1] にサポートを持ちます。次に、W_2n は、異なるサポート [0,1/2] および [1/2,1] を持つ、スケールが 1/2 のバージョンの 2 つの W_n を加算することで得られます。W_2n+1 は、同じバージョンの W_n を減算することで得られます。

n = 0 ～ 7 の場合、図Haar ウェーブレット パケットに示した W 関数が得られます。

Haar ウェーブレット パケット

これは次のコマンドを使用して取得できます。

[wfun,xgrid] = wpfun('db1',7,5);

ここで、wfun には、サポート xgrid の 1/2⁵ のグリッド上で計算した W_n (n = 0 ～ 7) の近似値が返されます。

元のウェーブレットとして、より正則なウェーブレットから始めて、同様の作成方法を使用すると、このような W 関数系の、より滑らかなバージョンが得られます。すべて、区間 [0, 2N–1] にサポートを持ちます。図db2 ウェーブレット パケットは、db2 ウェーブレットを元にした W 関数系を示しています。

db2 ウェーブレット パケット

ウェーブレット パケットの原子

関数 $$(W_n(x),n\in N)$$ から始めて、直交ウェーブレットを導くのと同じ方針に従い、3 つのインデックスの付いた解析関数のファミリ (波形) を検討します。

$$(W_{j,n,k}(x)=2^{-j/2}{W}_n(2^{-j}x-k)$$

ここで、n∊N、(j,k)∊Z² です。

ウェーブレットのフレームワークと同様に、k は時間の局所化パラメーターで、j はスケール パラメーターと解釈できます。では、n の解釈は何でしょうか。

ウェーブレット パケットの基本的な考え方は、j および k の固定値に対して、W_j,n,k はおおよそ 2^j· k の位置における信号の変動を、スケール 2^j で、最後のパラメーター n の許容される異なる値に対するさまざまな周波数で解析するというものです。

実際、Haar ウェーブレット パケットおよびdb2 ウェーブレット パケットに示されたウェーブレット パケットを注意深く確認すると、自然な順序の W_n (n = 0, 1, ..., 7) は、振動数で定義される順序とは厳密には一致しません。正確に言うと、db1 ウェーブレット パケットのゼロクロッシング (アップ クロッシングとダウン クロッシング) の回数を数えると、次のようになります。

自然順 n	0	1	2	3	4	5	6	7
db1 Wn のゼロクロッシングの回数	2	3	5	4	9	8	6	7

そのため、主な周波数がその順序により単調に増加するという性質を取り戻すため、自然順から再帰的に取得した "周波数順" を定義すると便利です。

自然順 n	0	1	2	3	4	5	6	7
周波数順 r(n)	0	1	3	2	6	7	5	4

前の図からわかるように、W_r(n)(x) は、おおよそ n 回 "振動" します。

信号 (たとえば、例 2 のチャープ) を解析するには、ウェーブレット パケット係数を自然順ではなく周波数順にしたがって、低い周波数 (一番下) から高い周波数 (一番上) へとプロットするのが適切です。

係数をプロットする場合、"周波数" 順か "自然" 順かの選択に関連するさまざまなオプションが、ウェーブレット アナライザー アプリを使用して利用できます。

これらのオプションは、関数 wpviewcf を使用すると、コマンド ライン モードからも利用できます。

ウェーブレット パケットの構成

関数の集合 W_j,n = (Wj,n,k(x),k∊Z) が (j,n) ウェーブレット パケットです。正の値の整数 j および n に対して、ウェーブレット パケットはツリー状に構成されます。図ツリー状に構成されたウェーブレット パケット (スケール j は深さを定義し、周波数 n はツリー内の位置を定義)のツリーは、最大レベルが 3 に等しい分解となるよう作成されています。各スケール j に対して、パラメーター n の取りうる値は 0,1, ...,2^j–1 です。

ツリー状に構成されたウェーブレット パケット (スケール j は深さを定義し、周波数 n はツリー内の位置を定義)

W_j,n (j はスケール パラメーター、n は周波数パラメーター) という表記は、通常のツリーの深さ-位置のラベル付けと一致しています。

ここで、$$W_{0,0}=(\varphi (x-k),k\in Z)$$、$$W_{1,1}=(\psi (\frac{x}{2}-k),k\in Z)$$ です。

ウェーブレット パケット基底のライブラリには、ウェーブレット基底に加え、その他の基底も複数含まれていることがわかります。そのような基底をいくつか見てみましょう。正確に説明するため、解析対象の信号が存在する (ファミリ W_0,0 が張る) 空間を V₀ とします。その場合、(W_d,1; d ≥ 1) は V₀ の直交基底です。

厳密に正である各整数 D に対して、(W_D,0, (W_d,1; 1 ≤ d ≤ D)) は V₀ の直交基底です。

また、関数のファミリ {(W_j+1,2n), (W_j+1,2n+1)} は、W_j,n が張る空間の直交基底であることもわかります。これは 2 つの部分空間に分割されます。W_j+1,2n が第 1 の部分空間を張り、W_j+1,2n+1 が第 2 の部分空間を張ります。

この最後の性質によって、ウェーブレット パケットの構成ツリーにおける分割を正確に解釈できます。作成されるノードはすべて、図ウェーブレット パケット ツリー: 分割と結合に示す形式だからです。

ウェーブレット パケット ツリー: 分割と結合

したがって、ツリー全体のうち、結合されている各二分部分木の葉が、最初の空間の直交基底に相当します。

V₀ に属する有限エネルギーの信号に対して、任意のウェーブレット パケット基底が、正確な再構成を提供します。また、周波数スケールのサブバンドでの情報割り当てを使用して、信号を符号化する固有の方法を提供します。

最適な分解の選択

ウェーブレット パケット ライブラリの構成に基づき、与えられた直交ウェーブレットから生じる分解の数を数えるのは自然なことです。

長さ N = 2^L の信号は、α 通りの方法で展開できます。ここで、α は深さ L の完全二分木の二分部分木の数です。その結果、$$\alpha \ge 2^{N/2}$$ となります ([Mal98] の 323 ページを参照)。

この数値は非常に大きい場合があり、また明示的な列挙は一般に取り扱いが困難であるため、効率的なアルゴリズムで計算できる便利な基準を使って最適な分解を見つけるというのは興味深いことです。最小限の基準については探求を続けています。

加法的性質を実証する関数は、二分木構造の効率的な検索と基本的な分割に適しています。古典的なエントロピーベースの基準はこれらの条件に適合し、与えられた信号の正確な表現について情報関連の性質を表します。エントロピーは、さまざまな分野に共通の概念であり、主に信号処理で使われます。ここでは 4 種類のエントロピー基準を紹介します ([CoiW92] を参照)。これ以外にも多数が利用でき、簡単に統合できます (「help wentropy」と入力)。以下の式で、s は信号、(s_i) は正規直交基底での s の係数です。

エントロピー E は、E(0) = 0 かつ次を満たす加法的コスト関数でなければなりません。

$$E(s)={\displaystyle {\sum }_{i}E(s_i)}$$

• (非正規化) Shannon エントロピー

  $$E1(s_i)=-s_i^2\log (s_i^2)$$

  したがって

  $$E1(s)=-{\displaystyle {\sum }_i{s}_i^2\log (s_i^2)}$$

  ただし、0log(0) = 0 とします。

• l^p ノルム (1 ℜ ≤ p) での集中度

  $$E2(s_i)={\left|s_i\right|}^p$$

  したがって

  $$E2(s)={\displaystyle {\sum }_i{\left|s_i\right|}^p}={\left|s\right|}_p^p$$

• "エネルギー" の対数エントロピー

  $$E3(s_i)=\log (s_i^2)$$

  したがって

  $$E3(s)={\displaystyle {\sum }_i\log (s_i^2)}$$

  ただし、log(0) = 0 とします。

• しきい値エントロピー

  $$\left|s_i\right|>\epsilon$$ では $$E4(s_i)=1$$ であり、それ以外では 0 であるため、E4(s) = # {i; $$\left|s_i\right|>\epsilon$$} は、信号がしきい値 ε より大きい時点の数です。

これらのエントロピー関数は、wentropy ファイルを使用して利用できます。

例 1: さまざまなエントロピーの計算

1. エネルギーが 1 に等しい信号を生成します。

   s = ones(1,16)*0.25;
   

2. s の Shannon エントロピーを計算します。

   e1 = wentropy(s,'shannon')
       e1 = 2.7726
   

3. s の l^1.5 エントロピーを計算します。これは norm(s,1.5)^1.5 と等価です。

   e2 = wentropy(s,'norm',1.5)
       e2 = 2
   

4. s の "対数エネルギー" エントロピーを計算します。

   e3 = wentropy(s,'log energy')
       e3 = -44.3614
   

5. s のしきい値エントロピーをしきい値 0.24 で計算します。

   e4 = wentropy(s,'threshold', 0.24)
       e4 = 16
   

例 2: 最小エントロピー分解

この簡単な例では、エントロピーを使用して、新しい分割が最小エントロピー分解を取得する上で関心の対象となるかどうかを判断する方法について説明します。

1. 元の信号は定数信号から始めます。2 つの情報があれば、この信号を定義および復元できます (長さと定数値)。

   w00 = ones(1,16)*0.25;
   

2. 元の信号のエントロピーを計算します。

   e00 = wentropy(w00,'shannon')
       e00 = 2.7726
   

3. 次に、Haar ウェーブレットを使用して w00 を分割します。

   [w10,w11] = dwt(w00,'db1');
   

4. レベル 1 の Approximation のエントロピーを計算します。

   e10 = wentropy(w10,'shannon')
       e10 = 2.0794
   

   レベル 1 の Detail w11 は 0 で、エントロピー e11 は 0 です。加法的性質のため、分解のエントロピーは e10+e11=2.0794 で与えられます。これを最初のエントロピー e00=2.7726 と比較しなければなりません。e10 + e11 < e00 なので、分割は関心の対象となります。

5. 今度は w10 を分割します (w11 でない理由は、エントロピーが 0 なのでヌルベクトルの分割は関心の対象外となるからです)。

   [w20,w21] = dwt(w10,'db1');
   

6. w20=0.5*ones(1,4) となり、w21 は 0 です。Approximation レベル 2 のエントロピーは次のようになります。

   e20 = wentropy(w20,'shannon')
       e20 = 1.3863
   

   再び e20 + 0 < e10 となるので、分割によってエントロピーは減少します。

7. 次になります。

   [w30,w31] = dwt(w20,'db1');
   e30 = wentropy(w30,'shannon')
       e30 = 0.6931
   
   [w40,w41] = dwt(w30,'db1')
       w40 = 1.0000
       w41 = 0
   
   e40 = wentropy(w40,'shannon')
       e40 = 0
   

   最後の分割処理において、1 つの情報のみが元の信号の再構成に必要であることがわかります。レベル 4 のウェーブレット基底は、Shannon エントロピーによると (e40+e41+e31+e21+e11 = 0 なのでエントロピーが最適な 0 となる) 最も良い基底です。

8. 例 1 で定義した信号 s のウェーブレット パケット分解を実行します。

   t = wpdec(s,4,'haar','shannon');
   

   エントロピー値のウェーブレット パケット ツリーは、ノードが元のエントロピーの数値でラベル付けされています。

   エントロピー値

   

9. ベスト ツリーを計算します。

   bt = besttree(t);
   

   ベスト ツリーを次の図に示します。この場合、ベスト ツリーはウェーブレット ツリーに相当します。ノードは最適エントロピーでラベル付けされています。

   最適エントロピー値

   

興味深い部分木

ウェーブレット パケットを使用するには、ツリーに関連する操作とラベル付けが必要です。ユーザー インターフェイスの実装は、これを考慮して構築されています。技術的な詳細については、リファレンス ページを参照してください。

レベル D で作成されたウェーブレット パケット分解ツリーに対応する深さ D の完全二分木は、WPT で示されます。

興味深い部分木として次のものがあります。

分解ツリー	葉の集合が基底となる部分木
ウェーブレット パケット分解ツリー	完全二分木: 深さ D の WPT
ウェーブレット パケット最適分解ツリー	WPT の二分部分木
ウェーブレット パケット ベスト レベル ツリー	WPT の完全二分部分木
ウェーブレット分解ツリー	深さ D の WPT の左片側二分部分木
ウェーブレット ベスト基底ツリー	WPT の左片側二分部分木

エントロピー基準 E に関して最適な分解の定義を以下のように導出しています。

分解	最適な分解	ベストレベル分解
ウェーブレット パケット分解	2D 個のツリーから検索	D 個のツリーから検索
ウェーブレット分解	D 個のツリーから検索	D 個のツリーから検索

終端以外のノードについて、次の基本手順を使用して、与えられたエントロピー基準 E に関して最適な部分木を見つけます (ここで、Eopt は最適なエントロピー値を示します)。

エントロピーの条件	ツリーおよびエントロピーのラベル付けに関する操作
$$E(node)\le {\displaystyle \sum _{c\text{ child of node}}Eopt(c)}$$	node≠root であれば、結合して、Eopt(node) = E(node) と設定
$$E(node)>{\displaystyle \sum _{c\text{ child of node}}Eopt(c)}$$	分割して、$$Eopt(node)={\displaystyle \sum _{c\text{ child of node}}Eopt(c)}$$ と設定

ただし、参照ツリーの固有の初期条件として、各終端ノード t に対して Eopt(t) = E(t) とします。

ノードからの信号 Approximation の再構成

関数 wprcoef を使用して、ウェーブレット パケット ツリーの任意のノードから信号の Approximation を再構成できます。これは、作業にウェーブレット パケット ツリー全体を使用しているか、最適性の基準で決定した部分木を使用しているかにかかわらず、当てはまります。信号の Approximation を再構成せずにノードからウェーブレット パケット係数を抽出する場合は、wpcoef を使用します。

ノイズがある Doppler 信号を読み込みます。

load noisdopp

sym4 ウェーブレットを使用して、レベル 5 までウェーブレット パケット分解を計算します。周期化モードを使用します。

dwtmode('per');
T = wpdec(noisdopp,5,'sym4');
plot(T)

ウェーブレット パケットの二分木をプロットし、(4,1) の組 (ノード 16) をクリックします。

ノード 16 からウェーブレット パケット係数を抽出します。

wpc = wpcoef(T,16);
% wpc is length 64

ノード 16 から信号の Approximation を取得します。

rwpc = wprcoef(T,16);
% rwpc is length 1024
plot(noisdopp,'k'); hold on;
plot(rwpc,'b','linewidth',2);
axis tight;

最適なウェーブレット パケットの二分木を確認します。

Topt = besttree(T);
% plot the best tree
plot(Topt)

(3,0) の組 (ノード 7) から信号の Approximation を再構成します。

rsig = wprcoef(Topt,7);
% rsig is length 1024
plot(noisdopp,'k'); hold on;
plot(rsig,'b','linewidth',2);
axis tight;

ウェーブレット パケットの二分木の中で抽出したい組がわかっている場合、depo2ind を使用してその組に対応するノードを特定できます。

たとえば、(3,0) の組に対応するノードを特定するには、次を入力します。

Node = depo2ind(2,[3 0]);

2 次元ウェーブレット パケット分解構造

ウェーブレット分解の場合とまったく同様に、前述の 1 次元フレームワークはイメージ解析に拡張できます。わずかな変更を直接行うと、四分木に関する定義になります。深さ 2 の例を次の図に示します。

深さ 2 の四分木

ウェーブレット パケットによる圧縮とノイズ除去

ウェーブレット パケットのフレームワークにおいて、圧縮とノイズ除去の考え方は、ウェーブレットのフレームワークで開発されたものと同じです。唯一新しい特徴は、より徹底的な解析によって柔軟性が向上することです。ウェーブレット パケットを使用した分解は 1 回で多数の基底を生成します。そこで、besttree とエントロピー関数を使用して、設計目的に対して最適な表現を探すことができます。