ewt

経験的ウェーブレット変換

R2020b 以降

構文

mra = ewt(x)

[mra,cfs] = ewt(x)

[mra,cfs,wfb] = ewt(x)

[mra,cfs,wfb,info] = ewt(x)

[___] = ewt(___,Name,Value)

ewt(___)

説明

mra = ewt(x) は、x の経験的ウェーブレット変換 (EWT) に対応する多重解像度解析 (MRA) 成分を返します。ewt は、経験的ウェーブレットとスケーリング フィルターを自動的に特定してエネルギーを維持する適応ウェーブレット細分化スキームを使用して信号を分解するために使用します。

既定では、経験的ウェーブレット フィルターの数は、x のマルチテーパー パワー スペクトル推定のピークを特定することによって自動的に決定されます。

[mra,cfs] = ewt(x) は、x の EWT 解析係数を返します。

[mra,cfs,wfb] = ewt(x) は、x の解析で使用される経験的ウェーブレット フィルター バンクを返します。

[mra,cfs,wfb,info] = ewt(x) は、x で特定されたピーク正規化周波数と、ウェーブレット フィルター バンクのおおよその周波数通過帯域を返します。

[___] = ewt(___,Name,Value) は、名前と値の組の引数を使用して追加オプションを指定します。これらの引数を前の入力構文のいずれかに追加できます。たとえば、'MaxNumPeaks',5 は、EWT フィルターの通過帯域を判定するために使用されるピークの最大数を 5 に指定します。

出力引数なしの ewt(___) は、経験的ウェーブレット MRA を使用して、元の信号を同じ Figure にプロットします。複素数値データの場合、実数部は MATLAB^® の色の順序の行列にある最初の色でプロットされ、虚数部は 2 番目の色でプロットされます。

例

経験的ウェーブレット変換の実行と信号のヒルベルト スペクトルの可視化

周波数にはっきりした変化が含まれる正弦波で構成される非定常連続信号を読み込み、可視化します。信号は 250 Hz でサンプリングされています。

fs = 250;
load nonstatdistinct
t = (0:length(nonstatdistinct)-1)/fs;
plot(t,nonstatdistinct)
xlabel('Time (s)')
ylabel('Signal')
axis tight

ewt を使用して、信号の多重解像度解析 (MRA) を取得します。

mra = ewt(nonstatdistinct);

関数hhtで MRA 成分を使用し、ヒルベルト スペクトルをプロットします。

hht(mra,fs)

周波数対時間のプロットは、MRA の各点における瞬間エネルギーを示す垂直方向のカラー バーがあるスパース プロットです。このプロットは、元の混合信号から分解された各成分の瞬時周波数スペクトルを表します。

経験的ウェーブレット変換フィルター バンクの可視化

周波数にはっきりした変化が含まれる正弦波で構成される非定常連続信号を作成します。信号は 1000 Hz でサンプリングされています。

Fs = 1000;
t = 0:1/Fs:4;
x1 = sin(2*pi*50*t) + sin(2*pi*200*t);
x2 = sin(2*pi*25*t) + sin(2*pi*100*t) + sin(2*pi*250*t);
x = [x1 x2] + 0.1*randn(1,length(t)*2);
t1 = (0:length(x)-1)/Fs;
plot(t1,x)
xlabel('Time (s)')
ylabel('Amplitude')
title('Signal')

ewt を使用して、信号の MRA を取得します。信号で特定されたピーク正規化周波数と、フィルター バンクのおおよその周波数通過帯域を表示します。周波数はサンプルあたりのサイクルであるため、サンプリング周波数で正規化します。ピーク周波数は正弦波の周波数に対応することに注意してください。

[mra,~,wfb,info] = ewt(x);
Fs*info.PeakFrequencies

ans = 5×1

  249.9375
  200.0750
  100.1000
   50.1125
   25.1187

Fs*info.FilterBank.Passbands 

ans = 5×2

  223.6941  500.0000
  141.5896  223.6941
   70.8573  141.5896
   35.4911   70.8573
         0   35.4911

信号の振幅スペクトルとフィルター バンクをプロットします。ピークの位置によってフィルターの通過帯域が決まります。

f = 0:Fs/length(x):Fs-1/length(x);
plot(f,wfb)
ylabel('Magnitude')
grid on
yyaxis right
plot(f,abs(fft(x)),'k--','linewidth',1.5)
ylabel('Magnitude')
xlabel('Hz')

経験的ウェーブレットはパーセバルのタイト フレームを形成するため、解析フィルター バンクと合成フィルター バンクは同じになります。したがって、各フィルターにおいて、各周波数の振幅を 2 乗したものの合計は 1 に等しくなります。合計が 1 に等しくない場合、完全再構成はできません。

ECG 信号の経験的ウェーブレット変換

ECG 信号を読み込みます。信号は 180 Hz でサンプリングされています。

load wecg

ewt を使用して、信号の多重解像度解析 (MRA) と対応する解析係数を取得します。4 つの最大ピークを使用して、フィルターの通過帯域を判定します。

mp = 4;
[mra,cfs] = ewt(wecg,'MaxNumPeaks',mp);

信号と MRA 成分をプロットします。

fs = 180;
subplot(mp+1,1,1)
t = (0:length(wecg)-1)/fs;
plot(t,wecg)
title('MRA of Signal')
ylabel('Signal')
axis tight
for k=1:mp
    subplot(mp+1,1,k+1)
    plot(t,mra(:,k))
    ylabel(['MRA ',num2str(k)])
    axis tight
end
xlabel('Time (s)')

MRA 成分を合計して、信号が完全再構成されるか検証します。

max(abs(wecg-sum(mra,2)))

ans = 8.8818e-16

EWT 解析係数のエネルギー維持を検証します。

cfsenergy = sum(sum(abs(cfs).^2));
[cfsenergy norm(wecg,2)^2]

ans = 1×2

  298.2759  298.2759

入力引数

x — 入力データ

ベクトル | timetable

入力データ。実数値または複素数値のベクトル、または single の列ベクトルを含む単一変数 timetable として指定します。x は少なくとも 2 つのサンプルをもたなければなりません。

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

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、コンマを使用して名前と値をそれぞれ区切り、Name を引用符で囲みます。

例: ewt(x,'MaxNumPeaks',5,'SegmentMethod','localmin') は、5 つの最大ピークと隣接するピーク間の最初の局所的最小値を使用して x の MRA を取得します。

PeakThresholdPercent — 最大ピークのしきい値の割合

70 (既定値) | 区間 (0, 100) 内の実数

x のマルチテーパー パワー スペクトルでどのピークを保持するかを決定するために使用される最大ピークのしきい値の割合。区間 (0, 100) 内の実数として指定します。x のマルチテーパー パワー スペクトル推定の局所的最大値は、最大ピークを 1 とする範囲 [0,1] に収まるように正規化されます。最大ピークが PeakThresholdPercent より厳密に大きい値になるすべてのピークが保持されます。

データ型: single | double

SegmentMethod — セグメンテーション手法

'geomean' (既定値) | 'localmin'

EWT フィルターの通過帯域を判定するために使用されるセグメンテーション手法。次として指定します。

• 'geomean' — 隣接するピークの幾何平均

• 'localmin' — 隣接するピーク間の最初の局所的最小値

隣接するピーク間で局所的最小値が特定されない場合、関数は幾何平均を使用します。

MaxNumPeaks — ピークの最大数

正の整数

EWT フィルターの通過帯域を判定するために使用されるピークの最大数。ewt が検出したピークが MaxNumPeaks で指定された数よりも少なかった場合、利用可能なピークの最大数が使用されます。ピークが見つからなかった場合、ewt はレベル 1 の離散ウェーブレット変換 (DWT) フィルター バンクを使用します。

MaxNumPeaks と PeakThresholdPercent の両方を指定することはできません。

データ型: single | double

FrequencyResolution — 周波数分解能帯域幅

0.25 以下の実数

マルチテーパー パワー スペクトル推定の周波数分解能帯域幅。0.25 以下の実数として指定します。

FrequencyResolution の値によって、マルチテーパー パワー スペクトル推定で使用される正弦テーパーの数が決まります。正弦マルチテーパー パワー スペクトル推定の帯域幅は (K+1)/(N+1) です。ここで、K はテーパーの数、N は信号の長さです。FrequencyResolution の最小値は 2.5/N です。ここで、N は信号長の最大値 64 です。

データ型: single | double

LogSpectrum — 対数スペクトル

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

ピーク周波数を判定するために対数スペクトルを使用するかどうかを示す logical 値。LogSpectrum が true に設定されている場合、対数のマルチテーパー パワー スペクトルが使用されます。PeakThresholdPercent のセグメンテーション手法を使用しており、他のピークよりも振幅が大幅に大きい支配的なピーク周波数が存在する場合は、LogSpectrum を true に設定することを検討してください。

出力引数

mra — 多重解像度解析

行列 | timetable

多重解像度解析 (MRA)。行列または timetable として返されます。

• x がベクトルの場合、mra は、抽出された MRA 成分が各列に格納された行列になります。

  • 実数値の x の場合、MRA 成分は中心周波数の降順に並べられます。mra の最後の列は、ローパス スケーリング フィルターに対応します。

  • 複素数値の x の場合、MRA 成分はサンプルあたり −½ サイクル付近から始まり、ローパス スケーリング係数が得られるまで中心周波数が減少します。その後、周波数はサンプルあたり +½ サイクルに向かって増加します。

• x が timetable の場合、mra は、複数の single 変数 (各変数に MRA 成分が格納される) をもつ timetable になります。

経験的ウェーブレットおよびスケーリング フィルターの周波数範囲の説明については、info 構造体配列を参照してください。

x が 64 サンプル未満の場合、ewt は長さ 64 になるように x をゼロ パディングしたものを使用します。MRA 成分は元の長さに切り捨てられます。

cfs — EWT 解析係数

行列

EWT 解析係数。行列として返されます。入力データが実数値の場合、cfs は実数値の行列になります。それ以外の場合、cfs は複素数値の行列になります。cfs の各列には、対応する MRA 成分の EWT 解析係数が格納されます。解析係数の周波数帯域は、MRA 成分の順序付けと同じです。x のサンプル数が 64 未満の場合、cfs には x のゼロ パディング バージョンから取得された解析係数が含まれます。

データ型: single | double

wfb — 経験的ウェーブレット フィルター バンク

行列

経験的ウェーブレット フィルター バンク。行列として返されます。wfb のフィルターの中心周波数は、mra および cfs の順序と一致します。経験的ウェーブレットはパーセバルのタイト フレームを形成するため、解析フィルター バンクと合成フィルター バンクは同じになります。したがって、MRA 成分の合計は、信号の完全再構成になります。

データ型: single | double

info — フィルター バンク情報

構造体配列

フィルター バンク情報。以下のフィールドをもつ構造体として返されます。

• PeakFrequencies — x で識別されるサイクル/サンプルのピーク正規化周波数 (列ベクトル)。実数値の x の場合、周波数は区間 (0, ½) で降順の正の値になります。複素数値の x の場合、周波数は (−½, ½) から順序付けられます。PeakFrequencies が空の場合、ewt はピークを検出せず、既定のレベル 1 の離散ウェーブレット変換 (DWT) 細分割が使用されます。

• FilterBank — MRAComponent と Passbands の 2 つの変数をもつ table。MRAComponent は、mra 内の MRA 成分の列インデックスです。Passbands は L 行 2 列の行列で、L は MRA 成分の数です。Passbands の各行は、対応する EWT フィルターと MRA 成分のおおよその周波数通過帯域 (サイクル/サンプル単位) です。

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

参照

[1] Gilles, Jérôme. “Empirical Wavelet Transform.” IEEE Transactions on Signal Processing 61, no. 16 (August 2013): 3999–4010. https://doi.org/10.1109/TSP.2013.2265222.

[2] Gilles, Jérôme, Giang Tran, and Stanley Osher. “2D Empirical Transforms. Wavelets, Ridgelets, and Curvelets Revisited.” SIAM Journal on Imaging Sciences 7, no. 1 (January 2014): 157–86. https://doi.org/10.1137/130923774.

[3] Gilles, Jérôme, and Kathryn Heal. “A Parameterless Scale-Space Approach to Find Meaningful Modes in Histograms — Application to Image and Spectrum Segmentation.” International Journal of Wavelets, Multiresolution and Information Processing 12, no. 06 (November 2014): 1450044. https://doi.org/10.1142/S0219691314500441.

拡張機能

C/C++ コード生成

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

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

• timetable の入力データはサポートされていません。

バージョン履歴

R2020b で導入

参考

アプリ

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

関数

• emd | modwtmra | vmd | hht

トピック

• Empirical Wavelet Transform
• Practical Introduction to Multiresolution Analysis
• 時間-周波数ギャラリー