Main Content

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

vmd

説明

imf = vmd(x)x の変分的モード分解を返します。vmd を使用して、複雑な信号をヒルベルト スペクトル解析に必要な有限数の固有モード関数 (IMF) に分解して単純化します。

[imf,residual] = vmd(x)x の変分的モード分解に対応する残差信号 residual を返します。

[imf,residual,info] = vmd(x) は、診断の目的のために IMF の追加情報 info と残差信号を返します。

[___] = vmd(x,Name,Value) では、1 つ以上の Name,Value のペア引数で指定された追加オプションを使用して、変分的モード分解を実行します。

vmd(___) では、同じ Figure 内で元の信号、IMF、残差信号をサブプロットとしてプロットします。

すべて折りたたむ

4 kHz でサンプリングされた信号を作成します。これはデジタル電話のすべてのキーをダイヤルした場合と似ています。信号を MATLAB® の timetable として保存します。

fs = 4e3;
t = 0:1/fs:0.5-1/fs;

ver = [697 770 852 941];
hor = [1209 1336 1477];

tones = [];

for k = 1:length(ver)
    for l = 1:length(hor)
        tone = sum(sin(2*pi*[ver(k);hor(l)].*t))';
        tones = [tones;tone;zeros(size(tone))];
    end
end

% To hear, type soundsc(tones,fs)

S = timetable(seconds(0:length(tones)-1)'/fs,tones);

timetable の変分的モード分解をプロットします。

vmd(S)

周波数 2 Hz、10 Hz、および 30 Hz の 3 つの正弦波からなる多成分信号を生成します。正弦波は 1 kHz で 2 秒間サンプリングされます。分散 0.01² のホワイト ガウス ノイズに信号を組み込みます。

fs = 1e3;
t = 1:1/fs:2-1/fs;
x = cos(2*pi*2*t) + 2*cos(2*pi*10*t) + 4*cos(2*pi*30*t) + 0.01*randn(1,length(t));

ノイズ信号の IMF を計算し、それを 3 次元プロットで可視化します。

imf = vmd(x);
[p,q] = ndgrid(t,1:size(imf,2));
plot3(p,q,imf)
grid on
xlabel('Time Values')
ylabel('Mode Number')
zlabel('Mode Amplitude')

計算された IMF を使用して、多成分信号のヒルベルト スペクトルをプロットします。周波数範囲を [0, 40] Hz に制限します。

hht(imf,fs,'FrequencyLimits',[0,40])

二次トレンド、チャープ、および t = 0.5 で 2 つの一定周波数間に急峻な遷移をもつ余弦波で構成される、区分合成信号を生成します。

x(t)=6t2+cos(4πt+10πt2)+{cos(60πt),t0.5,cos(100πt-10π),t>0.5.

信号は 1 kHz で 1 秒間サンプリングされます。個々の成分と合成信号をプロットします。

fs = 1e3;
t = 0:1/fs:1-1/fs;

x = 6*t.^2 + cos(4*pi*t+10*pi*t.^2) + ...
    [cos(60*pi*(t(t<=0.5))) cos(100*pi*(t(t>0.5)-10*pi))];

tiledlayout('flow')
nexttile
plot(t,[zeros(1,length(t)/2) cos(100*pi*(t(length(t)/2+1:end))-10*pi)])
xlabel('Time (s)')
ylabel('Cosine')

nexttile
plot(t,[cos(60*pi*(t(1:length(t)/2))) zeros(1,length(t)/2)])
xlabel('Time (s)')
ylabel('Cosine')

nexttile
plot(t,cos(4*pi*t+10*pi*t.^2))
xlabel('Time (s)')
ylabel('Chirp')


nexttile
plot(t,6*t.^2)
xlabel('Time (s)')
ylabel('Quadratic trend')

nexttile(5,[1 2])
plot(t,x)
xlabel('Time (s)')
ylabel('Signal')

変分的モード分解を実行して、4 つの固有モード関数を計算します。信号がもつ 4 つの異なる成分が復元されます。

[imf,res] = vmd(x,'NumIMFs',4);

tiledlayout('flow')

for i = 1:4
    nexttile
    plot(t,imf(:,i))
    txt = ['IMF',num2str(i)];
    ylabel(txt)
    xlabel('Time (s)')
    grid on
end

モード関数と残差を加算することにより、信号を再構成します。元の信号と再構成後の信号をプロットして比較します。

sig = sum(imf,2) + res;

nexttile(5,[1 2])
plot(t,sig,'LineWidth',1.5)

hold on

plot(t,x,':','LineWidth',2)
xlabel('Time (s)')
ylabel('Signal')
hold off
legend('Reconstructed signal','Original signal', ...
       'Location','northwest')

元の信号と再構成後の信号の差のノルムを計算します。

norm(x-sig',Inf)
ans = 0

この例では、MIT-BIH Arrhythmia Database の信号をラベル付けします [3]。データベース内の信号は 360 Hz でサンプリングされています。

記録 200 に対応する MIT データベース信号を読み込み、その信号をプロットします。

load mit200
Fs = 360;
plot(tm,ecgsig)
ylabel('Time (s)')
xlabel('Signal')

ECG 信号には、心拍リズムおよび振動する低周波パターンが生み出すスパイクが含まれます。ECG の個々のスポークは、重要な高次元の高調波を作成します。

ウィンドウを適用された信号の 9 つの固有モード関数を計算します。IMF を可視化します。

[imf,residual] = vmd(ecgsig,'NumIMF',9);
t = tiledlayout(3,3,'TileSpacing','compact','Padding','compact');
for n = 1:9
    ax(n) = nexttile(t);
    plot(tm,imf(:,n)')
    xlim([tm(1) tm(end)])
    txt = ['IMF',num2str(n)];
    title(txt)
    xlabel('Time (s)')
end
title(t,'Variational Mode Decomposition')

最初のモードは最もノイズが多く、2 番目のモードは心拍の頻度で振動しています。最初と最後の VMD モード以外を合計すると、低周波ベースライン振動と高周波ノイズの大部分が破棄されます。それにより、クリーンな ECG 信号を構築します。

cleanECG = sum(imf(:,2:8),2);
figure
plot(tm,ecgsig,tm,cleanECG)
legend('Original ECG','Clean ECG')
ylabel('Time (s)')
xlabel('Signal')

入力引数

すべて折りたたむ

等間隔でサンプリングされた時間領域信号。ベクトルまたは timetable として指定します。x が timetable の場合、増加する有限の行時間を含んでいなければなりません。timetable には、有限の負荷値をもつ数値データ ベクトルを 1 つのみ含めなければなりません。

メモ

timetable が欠損している場合や時間点が重複している場合、欠損または重複する時間および非等間隔の時間をもつ timetable の整理のヒントを使用して修正できます。

名前と値のペアの引数

オプションの引数 Name,Value のコンマ区切りペアを指定します。Name は引数名で、Value は対応する値です。Name は引用符で囲まなければなりません。Name1,Value1,...,NameN,ValueN のように、複数の名前と値のペアの引数を、任意の順番で指定できます。

例: 'NumIMF',10

モード収束絶対許容誤差。'AbsoluteTolerance' と正の実数スカラーで構成されるコンマ区切りのペアとして指定します。AbsoluteTolerance は最適化の停止条件の 1 つです。つまり、IMF の収束に向けた平均二乗の絶対的改善が 2 つの連続する反復において AbsoluteTolerance 未満となった場合に、最適化が停止します。

モード収束相対許容誤差。'RelativeTolerance' と正の実数スカラーで構成されるコンマ区切りのペアとして指定します。RelativeTolerance は最適化の停止条件の 1 つです。つまり、IMF の収束に向けた平均の相対的改善が 2 つの連続する反復において RelativeTolerance 未満となった場合に、最適化が停止します。

メモ

最適化プロセスは、'AbsoluteTolerance''RelativeTolerance' が一緒に達成された場合に停止します。

最適化反復の最大数。'MaxIterations' と正のスカラー整数で構成されるコンマ区切りのペアとして指定します。MaxIterations は最適化の停止条件の 1 つです。つまり、反復数が MaxIterations より大きくなった場合に、最適化が停止します。

MaxIterations は 0 または正の整数のみを使用して指定できます。

抽出された IMF の数。'NumIMF' と正のスカラー整数で構成されるコンマ区切りのペアとして指定します。

初期中心 IMF 周波数。'CentralFrequencies' と長さ NumIMFs のベクトルで構成されるコンマ区切りのペアとして指定します。ベクトル値は [0, 0.5] サイクル/サンプル内でなければなりません。これは、真の周波数が [0, fs/2] 内にあることを示します。ここで、fs はサンプル レートです。

初期 IMF。'InitialIMFs' および実数行列で構成されるコンマ区切りのペアとして指定します。行は時間サンプルに対応し、列はモードに対応します。

ペナルティ係数。'PenaltyFactor' と正の実数のスカラーで構成されるコンマ区切りのペアとして指定します。この引数により、再構成の忠実度が決まります。より厳しいデータ忠実度を得るには、ペナルティ係数の値を小さくします。

初期ラグランジュ乗数。'InitialLM' と複素数ベクトルで構成されるコンマ区切りのペアとして指定します。周波数領域における初期ラグランジュ乗数の範囲は、[0, 0.5] サイクル/サンプルです。この乗数は、再構成時に制約を設けます。乗数の長さは入力サイズに応じて異なります。

各反復におけるラグランジュ乗数の更新頻度。'LMUpdateRate' と正の実数スカラーで構成されるコンマ区切りのペアとして指定します。頻度が増せば収束は高速化しますが、最適化プロセスが局所的最適値で停止する可能性が高まります。

中心周波数の初期化方法。'InitializeMethod' と、'peaks''random' または 'grid' で構成されるコンマ区切りのペアとして指定します。

以下のように InitializeMethod を指定します。

  • 'peaks' は、周波数領域における信号のピーク位置として中心周波数を初期化します (既定)。

  • 'random' は、区間 [0,0.5] サイクル/サンプルで一様分布した乱数として中心周波数を初期化します。

  • 'grid' は、区間 [0,0.5] サイクル/サンプルで等間隔でサンプリングされたグリッドとして中心周波数を初期化します。

コマンド ウィンドウの進行状況表示の切り替え。'Display' と、'true' (もしくは 1) または 'false' (もしくは 0) で構成されるコンマ区切りのペアとして指定します。'true' を指定した場合、この関数は、モードの平均の絶対的および相対的な改善と 20 回の反復ごとの中心周波数を表示し、最後の停止情報を示します。

テーブルを表示するには Display に 1 を、非表示にするには 0 を指定します。

出力引数

すべて折りたたむ

固有モード関数。行列または timetable として返されます。各 imf は振幅変調および周波数変調された信号で、包絡線は正でゆっくりと変化します。各モードは、非減少で緩やかに変化し、中央値付近に集中する瞬時周波数をもちます。imf を使用してヒルベルト・ファン変換を適用し信号のスペクトル解析を実行します。

imf は以下として返されます。

  • x がベクトルの場合、各列が imf の行列

  • x が単一データ列の timetable の場合、timetable

残差信号。列ベクトルまたは単一データ列の timetable として返されます。residualvmd によって分解されていない元の信号 x の部分を表します。

residual は以下として返されます。

  • x がベクトルの場合、列ベクトル

  • x が単一データ列の timetable の場合、単一データ列の timetable

診断の詳細。次のフィールドがある構造体として返されます。

  • ExitFlag – 終了フラグ。値が 0 であれば、反復最大数に到達してアルゴリズムが停止したことを示します。値が 1 であれば、絶対許容誤差および相対許容誤差を満たしてアルゴリズムが停止したことを示します。

  • CentralFrequencies – IMF の中心周波数。

  • NumIterations – 合計反復回数。

  • AbsoluteImprovement – 最後の 2 回の反復間における IMF の収束に向けた平均二乗の絶対的改善。

  • RelativeImprovement – 最後の 2 回の反復間における IMF の収束に向けた平均の相対的改善。

  • LagrangeMultiplier – 最後の反復における周波数領域ラグランジュ乗数。

詳細

すべて折りたたむ

固有モード関数

関数 vmd は、信号 x(t)K 個という少数の狭帯域 "固有モード関数" (IMF) に分解します。

x(t)=k=1Kuk(t).

IMF は次の特性をもちます。

  1. 各モード uk は、次の形式をもつ振幅変調および周波数変調された信号です。

    uk(t)=Ak(t)cos(ϕk(t)),

    ここで、ϕk(t) はモードの位相であり、Ak(t) はその包絡線です。

  2. モードは、正でゆっくりと変化する包絡線をもちます。

  3. 各モードは、非減少で緩やかに変化し、中央値 fk 付近に集中する瞬時周波数 ϕ'k(t) をもちます。

変分的モード分解法では、すべてのモード波形とその中心周波数を同時に計算します。このプロセスは、制約のある変分の問題を最小化する uk(t)fk(t) のセットを見つけることで構成されています。

最適化

ukfk を計算するために、この手続きでは、次の最適な拡張ラグランジュ関数を見つけます。

L(uk(t),fk,λ(t))=αk=1Kddt[(δ(t)+jπt)uk(t)]ej2πfkt22+x(t)k=1Kuk(t)22+λ(t),x(t)k=1Kuk(t)(i)(ii)(iii),

ここで、内積 p(t),q(t)=p(t)q(t)dt かつ 2 ノルム p(t)22=p(t),p(t) です。正則化項 (i) には以下のステップが含まれます。

  1. ヒルベルト変換を使用して、各モードに関連する解析信号を計算します。ここで、* は畳み込みを表します。その結果、各モードは純粋に正のスペクトルをもちます。

  2. 解析信号を複素指数と乗算してベースバンドに復調します。

  3. 復調した解析信号がもつ勾配の 2 ノルムの二乗を計算することにより、帯域幅を推定します。

項 (ii) および (iii) は、二次ペナルティを適用し、ラグランジュ乗数を組み込むことにより、制約 x(t)=k=1Kuk(t) を適用します。PenaltyFactor α は、(ii) および (iii) と比較した (i) の相対的な重要度を測定します。

このアルゴリズムは、[1] で説明されている乗数の交互方向法を使用して、最適化問題を解決します。

アルゴリズム

関数 vmd は周波数領域の IMF を計算し、Uk(f) = DFT{uk(t)} により X(f) = DFT{x(t)} を再構成します。エッジの影響を除去するために、このアルゴリズムでは、信号のどちらかの側の長さを半分にミラーリングすることにより、その信号を拡張します。

最適化で紹介したラグランジュ乗数は、フーリエ変換 Ʌ(f) をもちます。ラグランジュ乗数ベクトルの長さは、拡張後の信号の長さになります。

'InitialIMFs' で特に指定されない限り、IMF はゼロに初期化されます。'InitializeMethod' で指定した方法のいずれかを使用して 'CentralFrequencies' を初期化します。vmd は、以下の条件のいずれかが達成されるまで、モードを反復的に更新します。

  • kukn+1(t)ukn(t)22/ukn(t)22<εr および kukn+1(t)ukn(t)22<εa が一緒に達成される。ここで、εr および εa は、それぞれ 'RelativeTolerance' および 'AbsoluteTolerance' を使用して指定されます。

  • アルゴリズムが、'MaxIterations' で指定した最大反復回数を超える。

(n + 1) 回目の反復の場合、アルゴリズムは次のステップを実行します。

  1. 次を計算するために、('NumIMF' を使用して指定された) 信号の K 個のモードを反復します。

    1. 次を使用して求めた、各モードの周波数領域波形。

      Ukn+1(f)=X(f)i<kUkn+1(f)i>kUkn(f)+Λn2(f)1+2α{2π(ffkn)}2,

      ここで、Ukn+1(f) は、(n + 1) 回目の反復で計算された k 番目のモードのフーリエ変換です。

    2. 次を使用して求めた、k 番目の中心周波数 fkn+1

      fkn+1=0|Ukn+1(f)|2fdf0|Ukn+1(f)|2dff|Ukn+1(f)|2|Ukn+1(f)|2.

  2. Λn+1(f)=Λn(f)+τ(X(f)kUkn+1(f)), を使用してラグランジュ乗数を更新します。ここで、τ は、'LMUpdateRate' を使用して指定される、ラグランジュ乗数の更新頻度です。

参照

[1] Boyd, Stephen, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. "Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers." Foundations and Trends® in Machine Learning. Vol 3, Number 1, 2011, pp. 1–122.

[2] Dragomiretskiy, Konstantin, and Dominique Zosso. "Variational Mode Decomposition." IEEE® Transactions on Signal Processing. Vol. 62, Number 3, 2014, pp. 531–534.

[3] Moody, George B., and Roger G. Mark. "The impact of the MIT-BIH Arrhythmia Database." IEEE Engineering in Medicine and Biology Magazine. Vol. 20, No. 3, May–June 2001, pp. 45–50.

参考

|

R2020a で導入