Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

modwt

最大重複離散ウェーブレット変換

説明

w = modwt(x) は、x の最大重複離散ウェーブレット変換 (MODWT) を返します。x は、実数値または複素数値のベクトルまたは行列です。x が行列の場合、modwtx の列を処理します。modwt がウェーブレット変換を計算するのは、x がベクトルの場合はレベル floor(log2(length(x))) まで、x が行列の場合はレベル floor(log2(size(x,1))) までです。既定では、modwt は、4 つの消失モーメント ('sym4') と周期的な境界の処理をもつ Daubechies 最小非対称ウェーブレットを使用します。

w = modwt(x,wname) は、MODWT に直交ウェーブレット wname を使用します。

w = modwt(x,Lo,Hi) は、スケーリング フィルター Lo とウェーブレット フィルター Hi を使用して、MODWT を計算します。これらのフィルターは、直交ウェーブレットの条件を満たさなければなりません。wname とフィルター ペア (LoHi) の両方を指定することはできません。

w = modwt(___,lev) は、前述の構文にある引数のいずれかを使用して、指定されたレベル lev まで MODWT を計算します。

w = modwt(___,'reflection') は、反射境界処理を使用して MODWT を計算します。他の入力は、前述の構文にある引数のいずれかです。ウェーブレット変換を計算する前に、modwt は、信号をその終端で信号長の 2 倍に対称的に拡張します。modwt が返すウェーブレット係数の数やスケーリング係数の数は、入力信号の長さの 2 倍に等しくなります。既定では、信号は周期的に拡張されます。

文字ベクトル 'reflection' 全体を入力しなければなりません。ウェーブレット マネジャーを使用して 'reflection' という名前のウェーブレットを追加した場合は、このオプションを使用する前に、そのウェーブレットの名前を変更しなければなりません。'reflection' は、入力引数リストで x より後の任意の位置に配置できます。

すべて折りたたむ

既定の sym4 ウェーブレットを最大レベルまで使用して、心電図 (ECG) 信号の MODWT を求めます。データは Percival & Walden (2000) の p.125 から採用しています (元々ワシントン大学の William Constantine と Per Reinhall によって提供されたデータ)。

load wecg;
wtecg = modwt(wecg);
whos wtecg
  Name        Size               Bytes  Class     Attributes

  wtecg      12x2048            196608  double              

wtecg の最初の 11 行は、スケール 21 から 211 のウェーブレット係数です。最後の行には、スケール 211 でのスケーリング係数が格納されます。スケール 23 の Detail (ウェーブレット) 係数をプロットします。

plot(wtecg(3,:))
title('Level 3 Wavelet Coefficients')

Figure contains an axes object. The axes object with title Level 3 Wavelet Coefficients contains an object of type line.

db2 ウェーブレットを最大レベルまで使用して、南方振動指数データの MODWT を求めます。

load soi;
wsoi = modwt(soi,'db2');

Fejér-Korovkin の長さ 8 のスケーリング フィルターとウェーブレット フィルターを使用して、ドイツ マルクとアメリカ ドルの為替レート データの MODWT を求めます。

load DM_USD;
[Lo,Hi] = wfilters('fk8');
wdm = modwt(DM_USD,Lo,Hi);

スケール 24 (レベル 4 に対応) までの ECG 信号の MODWT を求めます。既定の sym4 ウェーブレットを使用します。データは Percival & Walden (2000) の p.125 から採用しています (元々ワシントン大学の William Constantine と Per Reinhall によって提供されたデータ)。

load wecg;
wtecg = modwt(wecg,4);
whos wecg wtecg
  Name          Size              Bytes  Class     Attributes

  wecg       2048x1               16384  double              
  wtecg         5x2048            81920  double              

wtecg の行サイズは L+1 です。ここで、この場合、レベル (L) は 4 です。列サイズは入力サンプルの数と一致します。

反射境界処理を使用して、ECG 信号の MODWT を求めます。既定の sym4 ウェーブレットを使用して、レベル 4 までの変換を求めます。データは Percival & Walden (2000) の p.125 から採用しています (元々ワシントン大学の William Constantine と Per Reinhall によって提供されたデータ)。

load wecg;
wtecg = modwt(wecg,4,'reflection');
whos wecg wtecg
  Name          Size               Bytes  Class     Attributes

  wecg       2048x1                16384  double              
  wtecg         5x4096            163840  double              

wtecg は 4096 列で、これは入力信号の長さ wecg の 2 倍です。

23 チャネルの EEG データ Espiga3[3]を読み込みます。チャネルは列方向に配置されます。データは 200 Hz でサンプリングされます。

load Espiga3

最大重複離散ウェーブレット変換を最大レベルまで計算します。

wt = modwt(Espiga3);

信号エネルギーの 2 乗を求め、すべてのレベルのウェーブレット係数を合計して得られたエネルギーの 2 乗と比較します。ある成分のエネルギーが不釣り合いなほどに大きいため、エネルギーの 2 乗の対数を使用します。

sigN2 = vecnorm(Espiga3).^2;
wtN2 = sum(squeeze(vecnorm(wt,2,2).^2));
bar(1:23,log(sigN2))
hold on
scatter(1:23,log(wtN2),'filled','SizeData',100)
alpha(0.75)
legend('Signal Energy','Energy in Wavelet Coefficients', ...
        'Location','NorthWest')
xlabel('Channel')
ylabel('ln(squared energy)')

Figure contains an axes object. The axes object contains 2 objects of type bar, scatter. These objects represent Signal Energy, Energy in Wavelet Coefficients.

この例では、関数 MODWT と MODWTMRA の違いを示します。MODWT は、Detail 係数とスケーリング係数全体に信号のエネルギーを分割します。MODWTMRA はウェーブレット部分空間とスケーリング部分空間上に信号を投影します。

sym6 ウェーブレットを選択します。心電図 (ECG) 信号を読み込み、プロットします。ECG 信号のサンプリング周波数は 180 Hz です。データは Percival and Walden (2000) の p.125 から採用しています (元々 William Constantine and Per Reinhall, University of Washington によって提供されたデータ)。

load wecg
t = (0:numel(wecg)-1)/180;
wv = 'sym6';
plot(t,wecg)
grid on
title(['Signal Length = ',num2str(numel(wecg))])
xlabel('Time (s)')
ylabel('Amplitude')

Figure contains an axes object. The axes object with title Signal Length = 2048 contains an object of type line.

信号の MODWT を実行します。

wtecg = modwt(wecg,wv);

入力データは、N 個の時間点で評価された関数 f(x) のサンプルです。この関数は、さまざまなスケールと平行移動におけるスケーリング関数 ϕ(x) とウェーブレット ψ(x) の線形結合として f(x)=k=0N-1ck2-J0/2ϕ(2-J0x-k)+j=1J0fj(x) で表すことができます。ここで、fj(x)=k=0N-1dj,k2-j/2ψ(2-jx-k)J0 はウェーブレット分解のレベル数です。最初の合計は、粗いスケールによる信号の Approximation です。fj(x) は連続したスケールによる Detail です。MODWT は、展開の N 個の係数 {ck}(J0×N) 個の Detail 係数 {dj,k} を返します。wtecg の各行には、さまざまなスケールでの係数が含まれます。

長さ N の信号の MODWT を実行する場合、floor(log2(N)) レベルの分解があります (既定)。Detail 係数は各レベルで生成されます。スケーリング係数は最終レベルに対してのみ返されます。この例では、N=2048 なので、J0=floor(log2(2048))=11 になり、wtecg の行数は J0+1=11+1=12 になります。

MODWT は、||X||2=j=1J0||Wj||2+||VJ0||2 により、さまざまなスケールおよびスケーリング係数でエネルギーを分割します。ここで、X は入力データ、Wj はスケール j における Detail 係数、VJ0 は最終レベルのスケーリング係数です。

各スケールでエネルギーを計算し、その合計を評価します。

energy_by_scales = sum(wtecg.^2,2);
Levels = {'D1';'D2';'D3';'D4';'D5';'D6';...
    'D7';'D8';'D9';'D10';'D11';'A11'};
energy_table = table(Levels,energy_by_scales);
disp(energy_table)
    Levels     energy_by_scales
    _______    ________________

    {'D1' }         14.063     
    {'D2' }         20.612     
    {'D3' }         37.716     
    {'D4' }         25.123     
    {'D5' }         17.437     
    {'D6' }         8.9852     
    {'D7' }         1.2906     
    {'D8' }         4.7278     
    {'D9' }         12.205     
    {'D10'}         76.428     
    {'D11'}         76.268     
    {'A11'}         3.4192     
energy_total = varfun(@sum,energy_table(:,2))
energy_total=table
    sum_energy_by_scales
    ____________________

           298.28       

信号のエネルギーを計算し、それをすべてのスケールの合計エネルギーと比較することで、MODWT でエネルギーが維持されることを確認します。

energy_ecg = sum(wecg.^2);
max(abs(energy_total.sum_energy_by_scales-energy_ecg))
ans = 7.4402e-10

信号の MODWTMRA を実行します。

mraecg = modwtmra(wtecg,wv);

MODWTMRA は、さまざまなウェーブレット部分空間と最終的なスケーリング空間に対する関数 f(x) の投影を返します。つまり、MODWTMRA は、N 個の時間点で評価された k=0N-1ck2-J0/2ϕ(2-J0x-k)J0 個の {fj(x)} を返します。mraecg の各行は、f(x) を異なる部分空間に投影したものです。つまり、すべての投影を加算することで元の信号を復元することができます。これは、MODWT の場合は該当 "しません"wtecg の係数を加算しても、元の信号は復元 "されません"

時間点を選択し、その時間点で評価された f(x) の投影を加算し、元の信号と比較します。

time_point = 1000;
abs(sum(mraecg(:,time_point))-wecg(time_point))
ans = 3.0846e-13

MODWT とは異なり、MODWTMRA はエネルギーが維持される変換ではないことを確認します。

energy_ecg = sum(wecg.^2);
energy_mra_scales = sum(mraecg.^2,2);
energy_mra = sum(energy_mra_scales);
max(abs(energy_mra-energy_ecg))
ans = 115.7053

MODWTMRA は信号のゼロ位相フィルター処理です。特徴は時間で揃えられます。元の信号と投影の 1 つをプロットすることでこれを示します。配置をより適切に説明するには拡大します。

plot(t,wecg,'b')
hold on
plot(t,mraecg(4,:),'-')
hold off
grid on
xlim([4 8])
legend('Signal','Projection','Location','northwest')
xlabel('Time (s)')
ylabel('Amplitude')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Signal, Projection.

同じスケールで MODWT 係数を使用して同様のプロットを作成します。特徴は時間で揃えられないことに注意してください。MODWT は入力のゼロ位相フィルター処理では "ありません"

plot(t,wecg,'b')
hold on
plot(t,wtecg(4,:),'-')
hold off
grid on
xlim([4 8])
legend('Signal','Coefficients','Location','northwest')
xlabel('Time (s)')
ylabel('Amplitude')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Signal, Coefficients.

入力引数

すべて折りたたむ

ベクトルまたは行列として指定される入力信号。x がベクトルの場合、x には少なくとも 2 つの成分がなければなりません。x が行列の場合、x の行の次元は少なくとも 2 でなければなりません。

データ型: single | double
複素数のサポート: あり

解析ウェーブレット。次のいずれかとして指定します。

  • 'haar' — Haar ウェーブレット

  • 'dbN'N 個の消失モーメントをもつ極値位相 Daubechies ウェーブレット。ここで、N は 1 から 45 までの正の整数です。

  • 'symN'N 個の消失モーメントをもつ Symlet ウェーブレット。ここで、N は 2 から 45 までの正の整数です。

  • 'coifN'N 個の消失モーメントをもつ Coiflet ウェーブレット。ここで、N は 1 から 5 までの正の整数です。

  • 'fkN'N 個の係数をもつ Fejér-Korovkin ウェーブレット。ここで、N = 4, 6, 8, 14, 18, 22 です。

フィルター。偶数長の実数値ベクトルのペアとして指定します。Lo はスケーリング フィルター、Hi はウェーブレット フィルターです。フィルターは、直交ウェーブレットの条件を満たさなければなりません。LoHi の長さは等しくなければなりません。詳細については、wfilters を参照してください。ウェーブレット wname とフィルター ペア Lo,Hi の両方を指定することはできません。

変換レベル。floor(log2(N)) 以下の正の整数として指定します。ここで、x がベクトルの場合は N = length(x)x が行列の場合は N = size(x,1) です。lev を指定しない場合、既定で floor(log2(N)) になります。

出力引数

すべて折りたたむ

x の MODWT 変換。w には、x のウェーブレット係数と最終レベルのスケーリング係数が格納されます。x がベクトルの場合、wlev+1 行 N 列の行列です。x が行列の場合、wlev+1×N×NC の配列です。ここで、NC は x の列数です。'reflection' 境界処理を指定しない限り、N は入力信号長に等しくなります。指定した場合、N は入力信号長の 2 倍になります。配列 w の k 番目の行には、スケール 2k (ウェーブレット スケール 2(k-1)) のウェーブレット係数が格納されます。最後の (lev+1) 番目の行には、スケール 2lev のスケーリング係数が格納されます。

アルゴリズム

MODWT の標準アルゴリズムは、時間領域で循環畳み込みを直接実装します。MODWT のこの実装は、フーリエ領域で循環畳み込みを実行します。レベル j のウェーブレット フィルター係数とスケーリング フィルター係数は、離散フーリエ変換 (DFT) の積の逆離散フーリエ変換を取得することにより計算されます。積の DFT は、信号の DFT と j 番目のレベルのウェーブレット フィルターまたはスケーリング フィルターの DFT です。

Hk と Gk は、それぞれ MODWT のウェーブレット フィルターとスケーリング フィルターに対する長さ N の DFT を表します。j はレベルを表し、N はサンプル サイズを表します。

j 番目のレベルのウェーブレット フィルターは、

1Nk=0N1Hj,kei2πnk/N

によって定義されます。ここで

Hj,k=H2j1kmodNm=0j2G2mkmodN

です。

j 番目のレベルのスケーリング フィルターは

1Nk=0N1Gj,kei2πnk/N

です。ここで

Gj,k=m=0j1G2mkmodN

です。

参照

[1] Percival, Donald B., and Andrew T. Walden. Wavelet Methods for Time Series Analysis. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge ; New York: Cambridge University Press, 2000.

[2] Percival, Donald B., and Harold O. Mofjeld. “Analysis of Subtidal Coastal Sea Level Fluctuations Using Wavelets.” Journal of the American Statistical Association 92, no. 439 (September 1997): 868–80. https://doi.org/10.1080/01621459.1997.10474042.

[3] Mesa, Hector. “Adapted Wavelets for Pattern Detection.” In Progress in Pattern Recognition, Image Analysis and Applications, edited by Alberto Sanfeliu and Manuel Lazo Cortés, 3773:933–44. Berlin, Heidelberg: Springer Berlin Heidelberg, 2005. https://doi.org/10.1007/11578079_96.

拡張機能

バージョン履歴

R2015b で導入