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

modwt

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

構文

w = modwt(x)

w = modwt(x,wname)

w = modwt(x,Lo,Hi)

w = modwt(___,lev)

w = modwt(___,'reflection')

w = modwt(___,TimeAlign=alignflag)

説明

w = modwt(x) は、x の最大重複離散ウェーブレット変換 (MODWT) を返します。x は、実数値または複素数値のベクトルまたは行列です。x が行列の場合、modwt は x の列を処理します。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 とフィルター ペア (Lo と Hi) の両方を指定することはできません。

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

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

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

w = modwt(___,TimeAlign=alignflag) は、すべてのレベル (スケール) でウェーブレット係数とスケーリング係数を循環的にシフトして、スケーリング フィルターとウェーブレット フィルターの遅延を補正します。他の入力は、前述の構文にある引数のいずれかです。

例

既定のウェーブレットを使用する MODWT

既定の 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 行は、スケール $$2^1$$ から $$2^{11}$$ のウェーブレット係数です。最後の行には、スケール $$2^{11}$$ でのスケーリング係数が格納されます。スケール $$2^3$$ の Detail (ウェーブレット) 係数をプロットします。

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

2 つの消失モーメントをもつ Daubechies の極値位相ウェーブレットを使用する MODWT

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

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

スケーリング フィルターおよびウェーブレット フィルターを使用する MODWT

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

load DM_USD
[~,~,Lo,Hi] = wfilters("fk8");
wdmf = modwt(DM_USD,Lo,Hi);

同じウェーブレットで 2 番目の MODWT を取得しますが、今回はウェーブレット名を指定します。

wdmn = modwt(DM_USD,"fk8");

分解が等しいことを確認します。

max(abs(wdmf(:)-wdmn(:)))

ans = 0

指定したレベルまでの MODWT

スケール $$2^4$$ (レベル 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 です。列サイズは入力サンプルの数と一致します。

反射境界のある MODWT

反射境界処理を使用して、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 倍です。

MODWT 解析の時間揃え

ユニット インパルス信号を作成します。

n = 128;
sig = zeros(1,n);
sig(n/2) = 1;
clf
plot(sig)
axis tight
title("Unit Impulse")

modwt の既定の設定を使用して、信号の MODWT をレベル 4 まで取得します。時間揃え係数を使用して、信号の 2 番目の MODWT をレベル 4 まで取得します。

lev = 4;
w = modwt(sig,lev);
wTimeAligned = modwt(sig,lev,TimeAlign=true);

既定の設定で取得した MODWT のウェーブレット係数とスケーリング係数をプロットします。MODWT で使用されるウェーブレット フィルターとスケーリング フィルターの位相応答は非ゼロであるため、スケールと共に遅延が大きくなります。

for k=1:lev+1
    subplot(lev+1,1,k)
    plot(w(k,:))
    axis tight
    if k==1
        title("MODWT Analysis with Delay")
    end
end

時間を揃えた係数のプロットと比較します。

for k=1:lev+1
    subplot(lev+1,1,k)
    plot(wTimeAligned(k,:))
    axis tight
    if k==1
        title("Time-Aligned MODWT Analysis")
    end
end

多重信号の MODWT

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)')

MODWT と MODWTMRA の比較

この例では、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')

信号の MODWT を実行します。

wtecg = modwt(wecg,wv);

入力データは、$$N$$ 個の時間点で評価された関数 $$f(x)$$ のサンプルです。この関数は、さまざまなスケールと平行移動におけるスケーリング関数 $$\varphi (x)$$ とウェーブレット $$\psi (x)$$ の線形結合として $$f(x)=\sum _{k=0}^{N-1}{c}_k{2}^{-J_0/2}\varphi (2^{-J_0}x-k)+\sum _{j=1}^{J_0}{f}_j(x),$$ で表すことができます。ここで、$$f_j(x)=\sum _{k=0}^{N-1}{d}_{j,k}{2}^{-j/2}\psi (2^{-j}x-k)$$ と $$J_0$$ はウェーブレット分解のレベル数です。最初の合計は、粗いスケールによる信号の Approximation です。$$f_j(x)$$ は連続したスケールによる Detail です。MODWT は、展開の $$N$$ 個の係数 $$\{c_k\}$$ と $$(J_0\times N)$$ 個の Detail 係数 $$\{d_{j,k}\}$$ を返します。wtecg の各行には、さまざまなスケールでの係数が含まれます。

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

MODWT は、$${||X||}^2=\sum _{j=1}^{J_0}{||W_j||}^2+{||V_{J_0}||}^2,$$ により、さまざまなスケールおよびスケーリング係数でエネルギーを分割します。ここで、$$X$$ は入力データ、$$W_j$$ はスケール $$j$$ における Detail 係数、$$V_{J_0}$$ は最終レベルのスケーリング係数です。

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

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$$ 個の時間点で評価された $$\sum _{k=0}^{N-1}{c}_k{2}^{-J_0/2}\varphi (2^{-J_0}x-k)$$ と $$J_0$$ 個の $$\{f_j(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')

同じスケールで 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')

入力引数

x — 入力信号

ベクトル | 行列

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

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

wname — ウェーブレット

"sym4" (既定値) | 文字ベクトル | string スカラー

ウェーブレット。文字ベクトルまたは string スカラーとして指定します。ウェーブレットは直交でなければなりません。直交ウェーブレットは、ウェーブレット マネージャー wavemngr のタイプ 1 ウェーブレットとして指定されます。

有効な組み込み直交ウェーブレット ファミリは次のとおりです。最適局在化 Daubechies ("bl")、Beylkin ("beyl")、Coiflet ("coif")、Daubechies ("db")、Fejér-Korovkin ("fk")、Haar ("haar")、Han 線形位相モーメント ("han")、Morris 最小帯域幅 ("mb")、Symlet ("sym")、および Vaidyanathan ("vaid")。

各ファミリに含まれるウェーブレットの一覧については、wfilters を参照してください。ウェーブレット ファミリの略称を指定して waveinfo を使用することもできます。たとえば、waveinfo("db") のようにします。wavemngr("type",wn) を使用して、ウェーブレット wn が直交 (1 が返される) かどうかを判定します。たとえば、wavemngr("type","db6") は 1 を返します。

Lo,Hi — フィルター

偶数長の実数値のベクトル

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

メモ

数値パッケージでの modwt の実装における通常の慣例と一致するものとして、wfilters によって返される解析フィルターと合成フィルターの役割は、modwt では逆になります。スケーリング フィルターおよびウェーブレット フィルターを使用する MODWTを参照してください。

データ型: single | double

lev — 変換レベル

正の整数

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

alignflag — 時間揃え係数の logical 値

false または 0 (既定値) | true または 1

MODWT が、スケーリング フィルターとウェーブレット フィルターの遅延を補正するために、すべてのレベル (スケール) でウェーブレット係数とスケーリング係数を循環的にシフトするかどうかを決定する時間揃え係数の logical 値。数値または logical の 1 (true) または 0 (false) として指定します。係数は、Hess-Nielsen と Wickerhauser [4]の "エネルギー中心" 法を使用してシフトされます。係数のシフトは、信号内の特徴とウェーブレット係数の時間を揃えたい場合に便利です。

関数 imodwt を使用して信号を再構成する場合、または関数 modwtmra を使用して多重解像度解析を取得する場合は、係数をシフトしないでください。この場合、逆解析または多重解像度解析を取得する際に時間揃えが行われます。

データ型: logical

出力引数

w — MODWT 変換

行列 | 3 次元配列

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

アルゴリズム

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

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

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

$$\frac{1}{N}{\displaystyle \sum _{k=0}^{N-1}{H}_{j,k}}{e}^{i2\pi nk/N}$$

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

$$H_{j,k}=H_{2^{j-1}k\text{mod}N}{\displaystyle \prod _{m=0}^{j-2}{G}_{2^{m}k\text{mod}N}}$$

です。

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

$$\frac{1}{N}{\displaystyle \sum _{k=0}^{N-1}{G}_{j,k}}{e}^{i2\pi nk/N}$$

です。ここで

$$G_{j,k}={\displaystyle \prod _{m=0}^{j-1}{G}_{2^{m}k\text{mod}N}}$$

です。

参照

[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.

[4] Hess-Nielsen, N., and M.V. Wickerhauser. “Wavelets and Time-Frequency Analysis.” Proceedings of the IEEE 84, no. 4 (April 1996): 523–40. https://doi.org/10.1109/5.488698.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

• 入力 wname は定数でなければなりません。

GPU コード生成

GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

使用上の注意および制限:

• 入力 wname は定数でなければなりません。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

この関数は、GPU 配列を完全にサポートします。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2015b で導入

参考

アプリ

• 信号多重解像度アナライザー

関数

• dlmodwt | imodwt | modwtmra | modwtcorr | modwtvar | modwtxcorr

オブジェクト

• modwtLayer

トピック

• Practical Introduction to Multiresolution Analysis
• 時間-周波数ギャラリー
• Wavelet Analysis of Financial Data

×

例を開く

この例の変更されたバージョンがあります。編集された方の例を開きますか?

いいえ、変更されたバージョンを上書きします。 はい