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

emd

経験的モード分解

構文

[imf,residual] = emd(X)

[imf,residual,info] = emd(X)

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

emd(___)

説明

[imf,residual] = emd(X) は、X の経験的モード分解に対応する固有モード関数 imf と残差信号 residual を返します。emd を使用して、複雑な信号をヒルベルト スペクトル解析に必要な有限数の固有モード関数に分解して単純化します。

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

[___] = emd(___,Name,Value) では、1 つ以上の Name,Value のペア引数で指定された追加オプションを使用して、emd を推定します。

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

例

経験的モード分解の実行と信号のヒルベルト スペクトルの可視化

周波数にはっきりした変化が含まれる正弦波で構成される非定常連続信号を読み込み、可視化します。ジャックハンマーの振動と花火の音は、非定常連続信号の例です。信号は fs のレートでサンプリングされています。

load('sinusoidalSignalExampleData.mat','X','fs')
t = (0:length(X)-1)/fs;
plot(t,X)
xlabel('Time(s)')

混合信号には異なる振幅と周波数値をもつ正弦波が含まれます。

ヒルベルト スペクトル プロットを作成するには、信号の固有モード関数 (IMF) が必要です。経験的モード分解を実行して、IMF と信号の残差を計算します。信号が滑らかではないため、'pchip' を内挿法として指定します。

[imf,residual,info] = emd(X,'Interpolation','pchip');

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.026352   |  SiftMaxRelativeTolerance
      2      |        2     |    0.0039573   |  SiftMaxRelativeTolerance
      3      |        1     |     0.024838   |  SiftMaxRelativeTolerance
      4      |        2     |      0.05929   |  SiftMaxRelativeTolerance
      5      |        2     |      0.11317   |  SiftMaxRelativeTolerance
      6      |        2     |      0.12599   |  SiftMaxRelativeTolerance
      7      |        2     |      0.13802   |  SiftMaxRelativeTolerance
      8      |        3     |      0.15937   |  SiftMaxRelativeTolerance
      9      |        2     |      0.15923   |  SiftMaxRelativeTolerance
The decomposition stopped because the number of extrema of the residual signal is less than 'MaxNumExtrema'.

コマンド ウィンドウ内に生成されたテーブルは、生成された各 IMF のシフト反復の数、相対許容誤差、およびシフト停止基準を示します。この情報は info にも含まれます。名前と値のペア 'Display',0 を追加してテーブルを非表示にできます。

経験的モード分解を使用して取得した imf 成分を使用してヒルベルト スペクトル プロットを作成します。

hht(imf,fs)

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

正弦波の固有モード関数のゼロクロッシングと極値

この三角恒等式は、同じ物理量信号の 2 つの異なるビューを示します。

$\frac{5}{2}\cos \mathrm{2}\pi {\mathit{f}}_1\mathit{t}+\frac{1}{4}(\cos 2\pi \left({\mathit{f}}_1+{\mathit{f}}_2\right)t+\cos 2\pi \left({\mathit{f}}_1-{\mathit{f}}_2\right)t)=(2+\cos {}^2\pi {\mathit{f}}_2\mathit{t})\cos 2\pi {\mathit{f}}_1\mathit{t}$.

2 つの正弦波 s および z を生成します。s は 3 つの正弦波の和で、z は変調された振幅を持つ単一の正弦波です。これらの差の無限大ノルムを計算して、2 つの信号が等しいことを確認します。

t = 0:1e-3:10;
omega1 = 2*pi*100;
omega2 = 2*pi*20;
s = 0.25*cos((omega1-omega2)*t) + 2.5*cos(omega1*t) + 0.25*cos((omega1+omega2)*t);
z = (2+cos(omega2/2*t).^2).*cos(omega1*t);

norm(s-z,Inf) 

ans = 3.2729e-13

正弦波をプロットし、2 秒目から始まる 1 秒区間を選択します。

plot(t,[s' z'])
xlim([2 3])
xlabel('Time (s)')
ylabel('Signal')

信号のスペクトログラムを取得します。スペクトログラムは、3 つの異なる正弦波成分を示しています。フーリエ解析では、信号を正弦波の重ね合わせと見なします。

pspectrum(s,1000,'spectrogram','TimeResolution',4)

emd を使用して、信号の固有モード関数 (IMF) と追加の診断情報を計算します。既定では、この関数は各 IMF のシフティング反復の数、相対許容誤差、およびシフティング停止基準を示すテーブルを出力します。経験的モード分解では、信号を z と見なします。

[imf,~,info] = emd(s);

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        1     |   1.8025e-06   |  SiftMaxRelativeTolerance
The decomposition stopped because the current energy ratio is greater than 'MaxEnergyRatio'.

ゼロクロッシングと局所的極値の数は、最大で 1 だけ異なります。これは、信号が IMF であるために必要な条件を満たしています。

info.NumZerocrossing - info.NumExtrema

ans = 1

IMF をプロットし、2 秒目から始まる 0.5 秒区間を選択します。emd は信号を振幅変調された信号として表示するため、IMF は AM 信号です。

plot(t,imf)
xlim([2 2.5])
xlabel('Time (s)')
ylabel('IMF')

振動信号の固有モード関数の計算

破損したベアリングの振動信号をシミュレートします。経験的モード分解を実行して、信号の IMF を可視化し、欠陥を調査します。

ピッチの直径が 12 cm のベアリングは 8 つの回転要素を持ちます。各回転要素の直径は 2 cm です。内輪が 1 秒あたり 25 回駆動される間、外輪は静止状態を保ちます。加速度計はベアリングの振動を 10 kHz でサンプリングします。

fs = 10000;
f0 = 25;
n = 8;
d = 0.02;
p = 0.12;

正常なベアリングの振動信号には、駆動点周波数の次数が複数含まれます。

t = 0:1/fs:10-1/fs;
yHealthy = [1 0.5 0.2 0.1 0.05]*sin(2*pi*f0*[1 2 3 4 5]'.*t)/5;

共振は、測定プロセス中にベアリングの振動で励起されます。

yHealthy = (1+1./(1+linspace(-10,10,length(yHealthy)).^4)).*yHealthy;

共振によりベアリングの外輪に欠陥が生じることで、摩耗が進行します。欠陥があると、ベアリングの外輪転動体通過周波数 (BPFO) で繰り返される一連の影響を引き起こします。

$$BPFO=\frac{1}{2}n{f}_0[1-\frac{d}{p}\cos \theta ],$$

ここで、$$f_0$$ は駆動レート、$$n$$ は回転要素の数、$$d$$ は回転要素の直径、$$p$$ はベアリングのピッチの直径、$$\theta$$ はベアリングの接触角です。接触角は 15° と仮定して BPFO を計算します。

ca = 15;
bpfo = n*f0/2*(1-d/p*cosd(ca));

関数pulstranを使用して、影響を 5 ミリ秒の正弦波の周期列としてモデル化します。3 kHz の各正弦波に、フラット トップ ウィンドウによってウィンドウが適用されます。べき乗則を使用して、ベアリング振動信号に進行する摩耗を導入します。

fImpact = 3000;
tImpact = 0:1/fs:5e-3-1/fs;
wImpact = flattopwin(length(tImpact))'/10;
xImpact = sin(2*pi*fImpact*tImpact).*wImpact;

tx = 0:1/bpfo:t(end);
tx = [tx; 1.3.^tx-2];

nWear = 49000;
nSamples = 100000;
yImpact = pulstran(t,tx',xImpact,fs)/5;
yImpact = [zeros(1,nWear) yImpact(1,(nWear+1):nSamples)];

正常な信号に影響を追加して BPFO 振動信号を生成します。信号をプロットし、5.0 秒目から始まる 0.3 秒区間を選択します。

yBPFO = yImpact + yHealthy;

xLimLeft = 5.0;
xLimRight = 5.3;
yMin = -0.6;
yMax = 0.6;

plot(t,yBPFO)

hold on
[limLeft,limRight] = meshgrid([xLimLeft xLimRight],[yMin yMax]);
plot(limLeft,limRight,'--')
hold off

選択した区間にズームインして、影響の効果を可視化します。

xlim([xLimLeft xLimRight])

ホワイト ガウス ノイズを信号に付加します。$$1/150^2$$ のノイズ分散を指定します。

rn = 150;
yGood = yHealthy + randn(size(yHealthy))/rn;
yBad = yBPFO + randn(size(yHealthy))/rn;

plot(t,yGood,t,yBad)
xlim([xLimLeft xLimRight])
legend('Healthy','Damaged')

emd を使用して、正常なベアリング信号の経験的モード分解を実行します。最初の 5 つの固有モード関数 (IMF) を計算します。既定では、この関数は各 IMF のシフティング反復の数、相対許容誤差、およびシフティング停止基準を示すテーブルを出力します。

imfGood = emd(yGood,'MaxNumIMF',5);

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        3     |     0.017132   |  SiftMaxRelativeTolerance
      2      |        3     |      0.12694   |  SiftMaxRelativeTolerance
      3      |        6     |      0.14582   |  SiftMaxRelativeTolerance
      4      |        1     |     0.011082   |  SiftMaxRelativeTolerance
      5      |        2     |      0.03463   |  SiftMaxRelativeTolerance
The decomposition stopped because 'MaxNumIMF' was reached.

emd を出力引数なしで使用し、最初の 3 つのモードと残差を可視化します。'Display' を 0 に設定して、テーブルを非表示にします。

emd(yGood,'MaxNumIMF',5,'Display',0)

欠陥のあるベアリング信号の IMF を計算および可視化します。最初の経験的モードでは、高周波数に影響が見られます。この高周波数モードでは、摩耗が進行するにつれてエネルギーが増加します。3 番目のモードは、振動信号の共振を示します。

imfBad = emd(yBad,'MaxNumIMF',5);

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.041274   |  SiftMaxRelativeTolerance
      2      |        3     |      0.16695   |  SiftMaxRelativeTolerance
      3      |        3     |      0.18428   |  SiftMaxRelativeTolerance
      4      |        1     |     0.037177   |  SiftMaxRelativeTolerance
      5      |        2     |     0.095861   |  SiftMaxRelativeTolerance
The decomposition stopped because 'MaxNumIMF' was reached.

emd(yBad,'MaxNumIMF',5,'Display',0)

解析の次の手順は、抽出された IMF のヒルベルト スペクトルを計算することです。詳細については、振動信号のヒルベルト スペクトルの計算 (Signal Processing Toolbox)の例を参照してください。

信号の残差および固有モード関数の可視化

周波数にはっきりした変化が含まれる正弦波で構成される非定常連続信号を読み込み、可視化します。ジャックハンマーの振動と花火の音は、非定常連続信号の例です。信号は fs のレートでサンプリングされています。

load('sinusoidalSignalExampleData.mat','X','fs')
t = (0:length(X)-1)/fs;
plot(t,X)
xlabel('Time(s)')

混合信号には異なる振幅と周波数値をもつ正弦波が含まれます。

経験的モード分解を実行して、固有モード関数と信号の残差をプロットします。信号が滑らかではないため、'pchip' を内挿法として指定します。

emd(X,'Interpolation','pchip')

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.026352   |  SiftMaxRelativeTolerance
      2      |        2     |    0.0039573   |  SiftMaxRelativeTolerance
      3      |        1     |     0.024838   |  SiftMaxRelativeTolerance
      4      |        2     |      0.05929   |  SiftMaxRelativeTolerance
      5      |        2     |      0.11317   |  SiftMaxRelativeTolerance
      6      |        2     |      0.12599   |  SiftMaxRelativeTolerance
      7      |        2     |      0.13802   |  SiftMaxRelativeTolerance
      8      |        3     |      0.15937   |  SiftMaxRelativeTolerance
      9      |        2     |      0.15923   |  SiftMaxRelativeTolerance
The decomposition stopped because the number of extrema of the residual signal is less than 'MaxNumExtrema'.

元の信号、最初の 3 つの IMF、残差信号を含む対話型プロットが生成されます。コマンド ウィンドウ内に生成されたテーブルは、生成された各 IMF のシフト反復の数、相対許容誤差、およびシフト停止基準を示します。Display を 0 と指定することによってテーブルを非表示にできます。

プロット内の空白を右クリックし、[IMF セレクター] ウィンドウを開きます。[IMF セレクター] を使用して、生成された IMF、元の信号、および残差を選択的に表示します。

リストから表示する IMF を選択します。プロットに元の信号と残差を表示するかどうかを選択します。

これで選択した IMF がプロットに表示されます。

プロットを使用して、残差と共に元の信号から分解された個々の成分を可視化します。残差は IMF の総数に対して計算され、[IMF セレクター] ウィンドウで選択された IMF によって変化はしないことに注意してください。

入力引数

X — 等間隔でサンプリングされた時間領域信号

ベクトル | 単一データ列 timetable

一様にサンプリングされた時間領域。ベクトルまたは単一データ列 timetable として指定します。

名前と値のペアの引数

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

例: 'MaxNumIMF',5

'SiftRelativeTolerance' — コーシー型収束基準

0.2 (既定値) | 正のスカラー

コーシー型収束基準。'SiftRelativeTolerance' と正のスカラーで構成されるコンマ区切りのペアとして指定します。SiftRelativeTolerance はシフティング停止基準の 1 つです。つまり、現在の相対許容誤差が SiftRelativeTolerance 未満の場合にシフティングが停止します。

'SiftMaxIterations' — シフティング反復の最大数

100 (既定値) | 正のスカラー整数

シフティング反復の最大数。'SiftMaxIterations' と正のスカラー整数で構成されるコンマ区切りのペアとして指定します。SiftMaxIterations はシフティング停止基準の 1 つです。つまり、現在の反復数が SiftMaxIterations より大きい場合にシフティングが停止します。

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

'MaxNumIMF' — 抽出された IMF の最大数

10 (既定値) | 正のスカラー整数

抽出された IMF の最大数。'MaxNumIMF' と正のスカラー整数で構成されるコンマ区切りのペアとして指定します。MaxNumIMF は分解停止基準の 1 つです。つまり、生成された IMF の数が MaxNumIMF に等しい場合に分解が停止します。

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

'MaxNumExtrema' — 残差信号の極値の最大数

1 (既定値) | 正のスカラー整数

残差信号の極値の最大数。'MaxNumExtrema' と正のスカラー整数で構成されるコンマ区切りのペアとして指定します。MaxNumExtrema は分解停止基準の 1 つです。つまり、極値の数が MaxNumExtrema 未満の場合に分解が停止します。

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

'MaxEnergyRatio' — 信号対残差エネルギー比

20 (既定値) | スカラー

信号対残差エネルギー比。'MaxEnergyRatio' とスカラーで構成されるコンマ区切りのペアとして指定します。MaxEnergyRatio はシフティングの開始における信号のエネルギーと平均包絡線エネルギーとの比です。MaxEnergyRatio は分解停止基準の 1 つです。つまり、現在のエネルギー比が MaxEnergyRatio より大きい場合に分解が停止します。

'Interpolation' — 包絡線構築のための内挿法

'spline' (既定値) | 'pchip'

包絡線構築のための内挿法。'Interpolation' と 'spline' または 'pchip' のいずれかで構成されるコンマ区切りのペアとして指定します。

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

• X が滑らかな信号の場合は、'spline'

• X が滑らかでない信号の場合は、'pchip'

'spline' 内挿法は 3 次スプラインを使用します。一方で 'pchip' は区分的 3 次エルミート内挿多項式法を使用します。

'Display' — コマンド ウィンドウの情報表示の切り替え

1 (既定値) | 0

コマンド ウィンドウの情報表示の切り替えです。'Display' と 1 または 0 で構成されるコンマ区切りのペアで指定します。コマンド ウィンドウ内に生成されたテーブルは、生成された各 IMF のシフト反復の数、相対許容誤差、およびシフト停止基準を示します。テーブルを表示するには Display に 1 を、非表示にするには 0 を指定します。

出力引数

imf — 固有モード関数

行列 | timetable

固有モード関数。行列または timetable として返されます。imf は、包絡線がゼロに関して対称で、極値とゼロクロッシングの数の違いが多くても 1 つの関数です。imf を使用してヒルベルト・ファン変換を適用し信号のスペクトル解析を実行します。

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

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

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

residual — 信号の残差

列ベクトル | 単一データ列 timetable

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

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

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

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

info — 診断の詳細

構造体

診断の詳細。以下のフィールドを含む構造体として返します。

• NumIMF - 信号から抽出された IMF の数

• NumExtrema - 各 IMF における極値の数

• NumZeroCrossing - 各 IMF のゼロクロッシングの数

• NumSifting - 各 IMF の実行されたシフティングの数

• MeanEnvelopeEnergy - 各 IMF 計算において取得された上側および下側の包絡線の平均のエネルギー

• RelativeTolerance - 各 IMF の得られた相対許容誤差

アルゴリズム

Empirical Mode Decomposition

emd は、シフティング プロセスを使用して、信号 X(t) を k の数の固有モード関数 (IMF) と残差 r_k(t) に分解します。[1]と[2]にリストされているシフティング プロセスの簡単な概要は以下のとおりです。

1. 信号 X(t) の局所的最大値と局所的最小値を検索して上側包絡線 s₊(t) と、下側包絡線 s_–(t) を構築します。

2. i 番目の反復の平均包絡線、m_k,i(t) を計算します。

   $$m_{k,i}\left(t\right)=\frac{1}{2}\left[s_{+}\left(t\right)+s_{-}\left(t\right)\right]$$

3. 最初の反復については c_k(t) = X(t) を使用して、残差信号から平均包絡線を引きます。

   $$c_k\left(t\right)=c_k\left(t\right)-m_{k,i}\left(t\right)$$

   c_k(t) が IMF の基準に一致しない場合は、手順 4 および 5 がスキップされます。この手順を c_k(t) の新しい値で手順 1 において再度繰り返します。

4. c_k(t) が IMF の基準に一致する場合、新しい残差が計算されます。残差信号を更新するには、k 番目の IMF を前の残差信号から減算します。

   $$r_k\left(t\right)=r_{k-1}\left(t\right)-c_k\left(t\right)$$

5. 次に新しい信号 r_k(t) として取得した残差を使用して手順 1 から開始し、c_k(t) を固有モード関数として保存します。

固有モード関数 N の場合、元の信号は次のように表されます。

$$X\left(t\right)={\displaystyle \sum _{i=1}^N{\text{IMF}}_i\left(t\right)+r_N\left(t\right)}$$

シフティング プロセスの詳細は、[1]および[2]を参照してください。

SiftRelativeTolerance

SiftRelativeTolerance は[4]で提案されているコーシー型の停止基準です。現在の相対許容誤差が SiftRelativeTolerance 未満の場合にシフティングが停止します。現在の相対許容誤差は、以下のように定義されます。

$$\text{Relative Tolerance}\triangleq \frac{{\Vert r_{\text{prev}}\left(t\right)-r_{\text{cur}}\left(t\right)\Vert }_2^2}{{\Vert r_{\text{prev}}\left(t\right)\Vert }_2^2}.$$

コーシー基準がゼロクロッシングと局所的極値の数を直接カウントしないため、分解によって返された IMF は固有モード関数の厳しい定義を満たさない可能性があります。この場合、既定値から SiftRelativeTolerance の値を減らすことを試すことができます。停止条件の詳細は、[4]を参照してください。参照では、経験的モード分解で厳密に定義された IMF を求める利点と欠点についても説明します。

MaxEnergyRatio

エネルギー比は、シフティングの開始における信号のエネルギーと平均包絡線エネルギーの比です[3]。現在のエネルギー比が MaxEnergyRatio より大きい場合に分解が停止します。k IMF では、EnergyRatio は

$$\text{Energy Ratio}\triangleq 10{\log }_{10}\left(\frac{{\Vert X\left(t\right)\Vert }_2}{{\Vert r_i\left(t\right)\Vert }_2}\right).$$

と定義されます。

参照

[1] Huang, Norden E., Zheng Shen, Steven R. Long, Manli C. Wu, Hsing H. Shih, Quanan Zheng, Nai-Chyuan Yen, Chi Chao Tung, and Henry H. Liu. "The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis." Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences. Vol. 454, 1998, pp. 903–995.

[2] Rilling, G., Patrick Flandrin, and Paulo Gonçalves. "On empirical mode decomposition and its algorithms." IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing 2003. NSIP-03. Grado, Italy. 8–11.

[3] Rato, R.T., Manuel Ortigueira, and Arnaldo Batista. "On the HHT, its problems, and some solutions." Mechanical Systems and Signal Processing Vol. 22, 2008, pp. 1374–1394.

[4] Wang, Gang, Xian-Yao Chen, Fang-Li Qiao, Zhaohua Wu, and Norden Huang. "On Intrinsic Mode Function." Advances in Adaptive Data Analysis. Vol. 2, Number 3, 2010, pp. 277–293.

拡張機能

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

参考

hht

R2018a で導入

×

例を開く

この例の変更されたバージョンがシステム上にあります。代わりにこのバージョンを開きますか?

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