ウェーブレット パケット:Detail の分解

離散ウェーブレット変換および離散ウェーブレット パケット変換

次の図に、1 次元の入力信号の DWT ツリーを示します。

図 1: 1 次元入力信号のレベル 3 まで下げた DWT ツリー。

$$\underset{}{\overset{\sim }{G}}(f)$$ はスケーリング (ローパス) 解析フィルターであり、$$\underset{}{\overset{\sim }{H}}(f)$$ はウェーブレット (ハイパス) 解析フィルターを表します。下部のラベルは、周波数軸 [0,1/2] のサブバンドへの分割を示します。

図は、DWT の後続のレベルがローパス (スケーリング) フィルターの出力のみを処理することを示します。各レベルで、DWT は信号をオクターブ バンドに分割します。大きくサンプリングされている DWT では、バンドパス フィルターの出力は各レベルで 2 でダウンサンプリングされます。非間引き離散ウェーブレット変換では、出力はダウンサンプリングされません。

DWT ツリーと完全なウェーブレット パケット ツリーを比較します。

図 2: レベル 3 まで下げた完全なウェーブレット パケット ツリー。

ウェーブレット パケット変換では、フィルター処理はウェーブレット係数または Detail 係数にも適用されます。結果として、ウェーブレット パケットは、入力信号のサブバンド フィルター処理を徐々に細かく等幅間隔にします。各レベル $$j$$ で、周波数軸 [0,1/2] は $$2^j$$ 個のサブバンドに分割されます。レベル $$j$$ のヘルツ単位のサブバンドは、ほぼ以下のようになります。

$$[\frac{nFs}{2^{j+1}},\frac{(n+1)Fs}{2^{j+1}})n=0,1,\dots 2^j-1$$ ここで $$\text{Fs}$$ はサンプリング周波数です。

このバンドパスの Approximation がどの程度優れているかは、フィルターがどの程度周波数局在型であるかによって異なります。'fk18' (18 個の係数) などの Fejér-Korovkin フィルターの場合、Approximation は非常に良好です。Haar ('haar') などのフィルターの場合、Approximation は正確ではありません。

大きくサンプリングされているウェーブレット パケット変換では、バンドパス フィルターの出力は 2 でダウンサンプリングされます。非間引きウェーブレット パケット変換では、出力はダウンサンプリングされません。

ウェーブレット パケットを使用した時間-周波数解析

ウェーブレット パケットは DWT より周波数軸をより細かい間隔に分割するため、ウェーブレット パケットは時間-周波数解析より優れています。例として、加法性ノイズを伴う周波数 150 Hz および 200 Hz の 2 つの断続的な正弦波について考えます。データは 1 kHz でサンプリングされます。大きくサンプリングされているウェーブレット パケット変換でよくある時間分解能の損失を回避するには、非間引き変換を使用します。

dt = 0.001;
t = 0:dt:1-dt;
x = ...
cos(2*pi*150*t).*(t>=0.2 & t<0.4)+sin(2*pi*200*t).*(t>0.6 & t<0.9);
y = x+0.05*randn(size(t));
[wpt,~,F] = modwpt(x,'TimeAlign',true);
contour(t,F.*(1/dt),abs(wpt).^2)
grid on
xlabel('Time (secs)')
ylabel('Hz')
title('Time-Frequency Analysis -- Undecimated Wavelet Packet Transform')

ウェーブレット パケット変換は 150 Hz 成分と 200 Hz 成分を分離できることに注意してください。これが DWT ではできないのは、150 Hz と 200 Hz が同じオクターブ バンド内にあるためです。レベル 4 の DWT のオクターブ バンド (Hz 単位) を次に示します。

これらの 2 つの成分は DWT によって時間局在化されていますが、周波数で分離されていないことを示します。

mra  = modwtmra(modwt(y,'fk18',4),'fk18');
J = 5:-1:1;
freq = 1/2*(1000./2.^J);
contour(t,freq,flipud(abs(mra).^2))
grid on
ylim([0 500])
xlabel('Time (secs)')
ylabel('Hz')
title('Time-Frequency Analysis -- Undecimated Wavelet Transform')

もちろん、連続ウェーブレット変換 (CWT) は、ウェーブレット パケット変換と比べてより高い分解能の時間-周波数解析を提供しますが、ウェーブレット パケットには直交変換であること (直交ウェーブレットを使用する場合) によるさらなるメリットがあります。直交変換は、次の節で示すように、信号のエネルギーを維持して係数間で分割することを指します。

ウェーブレット パケット変換でのエネルギーの維持

ウェーブレット パケット変換は、信号を各レベルにおいて等幅サブバンドにフィルター処理するだけでなく、信号のエネルギーをサブバンド間で分割します。これは、間引きウェーブレット パケット変換と非間引きウェーブレット パケット変換の両方で行われます。これを示すために、地震の地震波記録を含むデータセットを読み込みます。このデータは、下記の信号分類の例で使用される時系列に似ています。

load kobe
plot(kobe)
grid on
xlabel('Seconds')
title('Kobe Earthquake Data')
axis tight

データの間引きウェーブレット パケット変換と非間引きウェーブレット パケット変換の両方をレベル 3 まで下げて取得します。間引きウェーブレット パケット変換での一貫した結果を確保するために、境界拡張モードを 'periodic' に設定します。

wptreeDecimated = dwpt(kobe,'Level',3,'Boundary','periodic');
wptUndecimated = modwpt(kobe,3);

間引きおよび非間引きの両方のウェーブレット パケットのレベル 3 係数で全エネルギーを計算し、元の信号のエネルギーと比較します。

decimatedEnergy = sum(cell2mat(cellfun(@(x) sum(abs(x).^2),wptreeDecimated,'UniformOutput',false)))

decimatedEnergy = 2.0189e+11

undecimatedEnergy = sum(sum(abs(wptUndecimated).^2,2))

undecimatedEnergy = 2.0189e+11

signalEnergy = norm(kobe,2)^2

signalEnergy = 2.0189e+11

DWT は、この重要な性質をウェーブレット パケット変換と共有します。

wt = modwt(kobe,'fk18',3);
undecimatedWTEnergy = sum(sum(abs(wt).^2,2))

undecimatedWTEnergy = 2.0189e+11

各レベルでのウェーブレット パケット変換によって周波数軸が等幅間隔に分割され、信号エネルギーがこれらの間隔に分割されるため、各ウェーブレット パケット内の相対的なエネルギーを維持するだけで、多くの場合に分類に便利な特徴ベクトルを作成できます。次の例では、これを示しています。

ウェーブレット パケット分類 -- 地震か爆発か

地震波記録は、多くのソースから活動を拾います。この活動をその出典に基づいて分類できることが重要です。たとえば、地震および爆発は同様の時間領域シグネチャを示す場合がありますが、これらの 2 つの出来事を区別することは非常に重要です。データセット EarthQuakeExplosionData には、2048 サンプルをもつ 16 の記録が含まれています。最初の 8 つの記録 (列) は地震の地震波記録であり、最後の 8 つの記録 (列) は爆発の地震波記録です。比較のためにデータを読み込み、地震および爆発の記録をプロットします。データは Shumway and Stoffer [3]を補完する R package astsa Stoffer [4]から取得されます。この例での使用については、著者の許可を得ています。

比較のために各グループの時系列をプロットします。

load EarthQuakeExplosionData
subplot(2,1,1)
plot(EarthQuakeExplosionData(:,3))
xlabel('Time')
title('Earthquake Recording')
grid on
axis tight
subplot(2,1,2)
plot(EarthQuakeExplosionData(:,9))
xlabel('Time')
title('Explosion Recording')
grid on
axis tight

非間引きウェーブレット パケット変換で 'fk6' ウェーブレットを使用して各時系列をレベル 3 まで下げて分解することで、ウェーブレット パケット特徴ベクトルを作成します。結果としておよその幅が 1/16 サイクル/サンプルのサブバンドが 8 つ作成されます。各サブバンドの相対エネルギーを使用して、特徴ベクトルを作成します。例として、最初の地震記録のウェーブレット パケットの相対エネルギー特徴ベクトルを取得します。

[wpt,~,F,E,RE] = modwpt(EarthQuakeExplosionData(:,1),3,'fk6');

RE には 8 つのサブバンドのそれぞれの相対エネルギーが含まれています。RE のすべての要素を合計すると、1 と等しくなります。これは大幅なデータ削減であることに注意してください。長さが 2048 の時系列は 8 つの要素をもつベクトルに削減され、各要素はレベル 3 のウェーブレット パケット ノードの相対エネルギーを表します。補助関数 helperEarthQuakeExplosionClassifier は、'fk6' ウェーブレットを使用してレベル 3 の 16 個の記録のそれぞれの相対エネルギーを取得します。結果として 16 行 8 列の行列が作成され、これらの特徴ベクトルを k-means 分類器への入力として使用します。分類器はギャップ統計の基準を使用して、1 から 6 の間の特徴ベクトルに最適なクラスター数を特定し、各ベクトルを分類します。さらに、分類器は、それぞれの時系列の 'fk6' ウェーブレットおよびパワー スペクトルを使用して取得したレベル 3 の非間引きウェーブレット変換係数に対してまったく同じ分類を実行します。非間引きウェーブレット変換によって、16 行 4 列の行列 (3 つのウェーブレット サブバンドと 1 つのスケーリング サブバンド) が作成されます。パワー スペクトルによって、16 行 1025 列の行列が作成されます。helperEarthQuakeExplosionClassifier のソース コードの一覧を付録に示します。分類器を実行するには、Statistics and Machine Learning Toolbox™ および Signal Processing Toolbox™ がなければなりません。

Level = 3;
[WavPacketCluster,WtCluster,PspecCluster] = ...
helperEarthQuakeExplosionClassifier(EarthQuakeExplosionData,Level)

WavPacketCluster = 
  GapEvaluation with properties:

    NumObservations: 16
         InspectedK: [1 2 3 4 5 6]
    CriterionValues: [-0.0039 0.0779 0.1314 0.0862 0.0296 -0.0096]
           OptimalK: 2

WtCluster = 
  GapEvaluation with properties:

    NumObservations: 16
         InspectedK: [1 2 3 4 5 6]
    CriterionValues: [-0.0095 0.0474 0.0706 0.0260 -0.0280 -0.0346]
           OptimalK: 2

PspecCluster = 
  GapEvaluation with properties:

    NumObservations: 16
         InspectedK: [1 2 3 4 5 6]
    CriterionValues: [0.3804 0.3574 0.3466 0.2879 0.3256 0.3326]
           OptimalK: 1

WavPacketCluster は、非間引きウェーブレット パケットの特徴ベクトルのクラスタリング出力であり、WtCluster は、非間引き DWT の特徴ベクトルのクラスタリング出力であり、PspecCluster は、パワー スペクトルに基づくクラスタリングの出力です。ウェーブレット パケット分類によって 2 つのクラスター (OptimalK: 2) が最適な数として特定されました。ウェーブレット パケット クラスタリングの結果を確認します。

res = [ WavPacketCluster.OptimalY(1:8)' ; WavPacketCluster.OptimalY(9:end)']

res = 2×8

     1     2     2     2     2     2     2     2
     1     1     1     1     1     1     1     2

エラーが 2 つしか生成されていないことがわかります。8 つの地震記録のうち、7 つはまとめて分類され (グループ 2)、1 つは誤ってグループ 1 に属すると分類されています。同じことが 8 つの爆発記録にも当てはまり、7 つがまとめて分類され、1 つが誤って分類されています。レベル 3 の非間引き DWT およびパワー スペクトルに基づく分類により、1 つのクラスターが最適解として返されます。

4 と等しいレベルを使用して上記を繰り返す場合、非間引きウェーブレット変換およびウェーブレット パケット変換は同様に実行されます。両方で 2 つのクラスターが最適なものとして特定され、両方で最初および最後の時系列が別のグループに属するものとして誤分類されます。

まとめ

この例では、ウェーブレット パケット変換と離散ウェーブレット変換の違いについて説明しました。特に、ウェーブレット パケット ツリーもウェーブレット (Detail) 係数のフィルター処理を行いますが、ウェーブレット変換が繰り返すのはスケーリング (Approximation) 係数のみです。

ウェーブレット パケット変換は、優れた周波数の分解能を提供する一方で、DWT のエネルギーが維持される重要な性質を共有することについて説明しました。一部の用途では、この優れた周波数の分解能により、ウェーブレット パケット変換が DWT の魅力的な代替手段となります。

付録

helperEarthQuakeExplosionClassifier

function [WavPacketCluster,WtCluster,PspecCluster] = helperEarthQuakeExplosionClassifier(Data,Level)  
%   This function is provided solely in support of the featured example
%   waveletpacketsdemo.mlx. 
%   This function may be changed or removed in a future release.
%   Data - The data matrix with individual time series as column vectors.
%   Level - The level of the wavelet packet and wavelet transforms

NP = 2^Level;
NW = Level+1;
WpMatrix = zeros(16,NP);
WtMatrix = zeros(16,NW);

for jj = 1:8
    [~,~,~,~,RE] = modwpt(Data(:,jj),Level,'fk6');
    wt = modwt(Data(:,jj),Level,'fk6');
    WtMatrix(jj,:) = sum(abs(wt).^2,2)./norm(Data(:,jj),2)^2;
    WpMatrix(jj,:) = RE;
end

for kk = 9:16
    [~,~,~,~,RE] = modwpt(Data(:,kk),Level,'fk6');
    wt = modwt(Data(:,kk),Level,'fk6');
    WtMatrix(kk,:) = sum(abs(wt).^2,2)./norm(Data(:,kk),2)^2;
    WpMatrix(kk,:) = RE;
end

% For repeatability
rng('default');

WavPacketCluster = evalclusters(WpMatrix,'kmeans','gap','KList',1:6,...
    'Distance','cityblock');

WtCluster = evalclusters(WtMatrix,'kmeans','gap','KList',1:6,...
'Distance','cityblock');

Pxx = periodogram(Data,hamming(numel(Data(:,1))),[],1,'power');

PspecCluster = evalclusters(Pxx','kmeans','gap','KList',1:6,...
    'Distance','cityblock');
end