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

hht

ヒルベルト・ファン変換

構文

hs = hht(imf)

hs = hht(imf,fs)

[hs,f,t] = hht(___)

[hs,f,t,imfinsf,imfinse] = hht(___)

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

hht(___)

hht(___,freqlocation)

説明

hs = hht(imf) は、固有モード関数 imf によって指定される信号のヒルベルト スペクトル hs を返します。hs は、スペクトル成分が経時変化する信号の混成から成る信号の解析に役立ちます。hht を使用して信号のヒルベルトスペクトル解析を実行し、局所化された特徴を特定します。

hs = hht(imf,fs) は、fs のレートでサンプリングされた信号のヒルベルト スペクトル hs を返します。

[hs,f,t] = hht(___) は、hs に加えて周波数ベクトル f と時間ベクトル t を返します。これらの出力引数は、前のいずれの入力構文でも使用できます。

[hs,f,t,imfinsf,imfinse] = hht(___) は、信号の診断のために、固有モード関数の瞬時周波数 imfinsf と瞬間エネルギー imfinse も返します。

[___] = hht(___,Name,Value) では、1 つ以上の Name,Value のペア引数で指定された追加オプションを使用して、ヒルベルト スペクトル パラメーターを推定します。

出力引数なしで hht(___) を使用すると、現在の Figure ウィンドウにヒルベルト スペクトルがプロットされます。この構文は、これより前の構文の任意の入力引数で使用できます。

hht(___,freqlocation) では、周波数軸の場所を指定するためにオプションの freqlocation 引数を使用してヒルベルト スペクトルをプロットします。既定の設定では、周波数は y 軸上に表されます。

例

二次チャープのヒルベルト スペクトル

ガウス変調二次チャープを生成します。2 kHz のサンプルレートと 2 秒の信号持続時間を指定します。

fs = 2000;
t = 0:1/fs:2-1/fs;
q = chirp(t-2,4,1/2,6,'quadratic',100,'convex').*exp(-4*(t-1).^2);
plot(t,q)

emd を使用して、固有モード関数 (IMF) と残差を可視化します。

emd(q)

信号の IMF を計算します。'Display' の名前と値のペアを使用して、各 IMF のふるい分け反復の数、相対許容誤差、およびふるい分け停止基準を示すテーブルを出力します。

imf = emd(q,'Display',1);

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |    0.0063952   |  SiftMaxRelativeTolerance
      2      |        2     |       0.1007   |  SiftMaxRelativeTolerance
      3      |        2     |      0.01189   |  SiftMaxRelativeTolerance
      4      |        2     |    0.0075124   |  SiftMaxRelativeTolerance
Decomposition stopped because the number of extrema in the residual signal is less than the 'MaxNumExtrema' value.

計算された IMF を使用して、二次チャープのヒルベルト スペクトルをプロットします。周波数範囲を 0 Hz ～ 20 Hz に制限します。

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

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

周波数にはっきりした変化が含まれる正弦波で構成される非定常連続信号を読み込み、可視化します。ジャックハンマーの振動と花火の音は、非定常連続信号の例です。信号は 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');

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

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

hht(imf,fs)

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

クジラの歌のヒルベルト スペクトル

4 kHz でサンプリングされた太平洋のシロナガスクジラのオーディオ データを含むファイルを読み込みます。ファイルは、コーネル大学の生物音響学研究プログラムが管理する動物発声ライブラリのものです。データの時間スケールは、音の高さを上げ鳴き声を聞き取りやすくするために係数 10 で圧縮されています。信号を MATLAB® の timetable に変換し、プロットします。信号のノイズの中で 4 つの特徴が目立っています。最初は "ふるえ声"、他の 3 つは "うめき声" として知られています。

[w,fs] = audioread('bluewhale.wav');
whale = timetable(w,'SampleRate',fs);
stackedplot(whale);

emd を使用して、最初の 3 つの固有モード関数 (IMF) と残差を可視化します。

emd(whale,'MaxNumIMF',3)

信号の最初の 3 つの IMF を計算します。'Display' の名前と値のペアを使用して、各 IMF のふるい分け反復の数、相対許容誤差、およびふるい分け停止基準を示すテーブルを出力します。

imf = emd(whale,'MaxNumIMF',3,'Display',1);

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        1     |      0.13523   |  SiftMaxRelativeTolerance
      2      |        2     |     0.030198   |  SiftMaxRelativeTolerance
      3      |        2     |      0.01908   |  SiftMaxRelativeTolerance
Decomposition stopped because maximum number of intrinsic mode functions was extracted.

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

hht(imf,'FrequencyLimits',[0 1400])

同じ周波数の範囲のヒルベルト スペクトルを計算します。ふるえ声とうめき声のヒルベルト スペクトルをメッシュ プロットとして可視化します。

[hs,f,t] = hht(imf,'FrequencyLimits',[0 1400]);

mesh(seconds(t),f,hs,'EdgeColor','none','FaceColor','interp')
xlabel('Time (s)')
ylabel('Frequency (Hz)')
zlabel('Instantaneous Energy')

信号のヒルベルト スペクトル パラメーターの計算

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

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

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

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

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

コマンド ウィンドウ内に生成されたテーブルは、生成された各 IMF のふるい分け反復の数、相対許容誤差、およびふるい分け停止基準を示します。この情報は info にも含まれます。'Display' を 0 と指定することによってテーブルを非表示にできます。

ヒルベルト スペクトルの次のパラメーターを計算します。ヒルベルト スペクトル hs、周波数ベクトル f、時間ベクトル t、瞬時周波数 imfinsf および瞬間エネルギー imfinse。

[hs,f,t,imfinsf,imfinse] = hht(imf,fs);

計算したヒルベルト スペクトル パラメーターを時間-周波数解析と信号の診断に使用します。

多成分信号の VMD

周波数 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])

振動信号のヒルベルト スペクトルの計算

破損したベアリングの振動信号をシミュレートします。この信号のヒルベルト スペクトルを計算し、欠陥を探します。

ピッチの直径が 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) を計算します。'Display' の名前と値のペアを使用して、各 IMF のふるい分け反復の数、相対許容誤差、およびふるい分け停止基準を示すテーブルを出力します。

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

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        3     |     0.016207   |  SiftMaxRelativeTolerance
      2      |        3     |     0.081827   |  SiftMaxRelativeTolerance
      3      |        6     |      0.14991   |  SiftMaxRelativeTolerance
      4      |        1     |     0.011853   |  SiftMaxRelativeTolerance
      5      |        2     |     0.015302   |  SiftMaxRelativeTolerance
Decomposition stopped because maximum number of intrinsic mode functions was extracted.

emd を出力引数なしで使用し、最初の 3 つのモードと残差を可視化します。

emd(yGood,'MaxNumIMF',5)

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

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

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.042507   |  SiftMaxRelativeTolerance
      2      |        4     |     0.011889   |  SiftMaxRelativeTolerance
      3      |        4     |      0.19587   |  SiftMaxRelativeTolerance
      4      |        1     |     0.016427   |  SiftMaxRelativeTolerance
      5      |        2     |     0.026917   |  SiftMaxRelativeTolerance
Decomposition stopped because maximum number of intrinsic mode functions was extracted.

emd(yBad,'MaxNumIMF',5)

欠陥のあるベアリング信号について、最初の経験的モードのヒルベルト スペクトルをプロットします。最初のモードは、高周波数の影響の効果をとらえています。ベアリングの摩耗が進行するにつれて、影響のエネルギーが増加します。

figure
hht(imfBad(:,1),fs)

3 番目のモードのヒルベルト スペクトルは、振動信号の共振を示します。周波数範囲を 0 Hz ～ 100 Hz に制限します。

hht(imfBad(:,3),fs,'FrequencyLimits',[0 100])

比較のために、正常なベアリング信号について、最初と 3 番目のモードのヒルベルト スペクトルをプロットします。

subplot(2,1,1)
hht(imfGood(:,1),fs)
subplot(2,1,2)
hht(imfGood(:,3),fs,'FrequencyLimits',[0 100])

入力引数

imf — 固有モード関数

行列 | timetable

固有モード関数。行列または timetable として指定します。imf は、包絡線がゼロに関して対称であり、極値とゼロクロッシングの数の違いが多くても 1 つの信号です。emd は、複雑な信号をヒルベルト スペクトル解析を実行するために必要な固有モード関数の有限数に分解および単純化するために使用されます。

hht は、imf の各列を固有モード関数として処理します。imf の計算の詳細は、emd を参照してください。

fs — サンプルレート

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

サンプルレート。正のスカラーで指定します。fs が指定されない場合、2π の正規化された周波数がヒルベルト スペクトルの計算に使用されます。imf が timetable として指定されている場合、これからサンプルレートが推定されます。

freqlocation — プロット上の周波数軸の位置

'yaxis' (既定値) | 'xaxis'

プロット上の周波数軸の位置。'yaxis' または 'xaxis' として指定します。プロットの y 軸または x 軸上に周波数データを表示するには、それぞれ 'yaxis' または 'xaxis' として freqlocation を指定します。

名前と値のペアの引数

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

例: 'FrequencyResolution',1

'FrequencyLimits' — ヒルベルト スペクトルを計算する周波数限界

[0,fs/2] (既定値) | 1 行 2 列の整数値のベクトル

ヒルベルト スペクトルを計算する周波数限界。'FrequencyLimits' と 1 行 2 列の整数値ベクトルで構成されるコンマ区切りのペアとして指定します。FrequencyLimits は Hz で指定します。

'FrequencyResolution' — 周波数範囲を離散化する周波数分解能

(f_high-f_low)/100 (既定値) | 正のスカラー

周波数限界を離散化する周波数分解能。'FrequencyResolution' と正のスカラーとで構成されるコンマ区切りのペアとして指定します。

FrequencyResolution は Hz で指定します。'FrequencyResolution' が指定されない場合、(f_high-f_low)/100 の値は FrequencyLimits から推定されます。ここで、f_high は FrequencyLimits の上限、f_low は下限です。

'MinThreshold' — ヒルベルト スペクトルの最小しきい値

-inf (既定値) | スカラー

ヒルベルト スペクトルの最小しきい値。'MinThreshold' とスカラーで構成されるコンマ区切りのペアとして指定します。

MinThreshold は、$$10{\log }_{10}\left(hs\right)$$ の対応する要素が MinThreshold より小さい場合に hs の要素を 0 に設定します。

出力引数

hs — 信号のヒルベルト スペクトル

スパース行列

信号のヒルベルト スペクトル。スパース行列として返されます。hs を使用して時間-周波数を解析し、信号の局所的な特徴を特定します。

f — 周波数値

ベクトル

信号の周波数値。ベクトルとして返されます。hht は、周波数ベクトル f と時間ベクトル t を使用してヒルベルト スペクトル プロットを作成します。

数学的には、f は f = f_low : f_res : f_high と表されます。ここで、f_res は周波数分解能です。

t — 時間値

ベクトル | duration 配列

信号の時間値。ベクトルまたは duration 配列として返されます。hht は、時間ベクトル t と周波数ベクトル f を使用してヒルベルト スペクトル プロットを作成します。

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

• imf を配列として指定した場合は配列。

• imf を一様にサンプリングされた timetable として指定した場合は duration 配列。

imfinsf — 各 IMF の瞬時周波数

ベクトル | 行列 | timetable

各 IMF の瞬時周波数。ベクトル、行列、または timetable として返されます。

imfinsf は imf と同じ列数であり、以下として返されます。

• imf をベクトルとして指定した場合はベクトル。

• imf を行列として指定した場合は行列。

• imf を一様にサンプリングされた timetable として指定した場合は timetable。

imfinse — 各 IMF の瞬間エネルギー

ベクトル | 行列 | timetable

各 IMF の瞬間エネルギー。ベクトル、行列、または timetable として返されます。

imfinse は imf と同じ列数であり、以下として返されます。

• imf をベクトルとして指定した場合はベクトル。

• imf を行列として指定した場合は行列。

• imf を一様にサンプリングされた timetable として指定した場合は timetable。

アルゴリズム

ヒルベルト・ファン変換は、非定常および非線形データの時間-周波数解析の実行に役立ちます。ヒルベルト・ファン手順は以下のステップから構成されます。

1. emd はデータセット x を固有モード関数の有限個に分解します。

2. 各固有モード関数 x_i の場合、関数 hht は次のとおりです。

   1. hilbert を使用して解析信号 $$z_i(t)=x_i(t)+jH\{x_i(t)\}$$ を計算します。ここで、H{x_i} は x_i のヒルベルト変換です。

   2. z_i を $$z_i(t)=a_i(t)e^{j{\theta }_i(t)}$$ として表します。ここで、a_i(t) は瞬間的な振幅で、$${\theta }_i(t)$$ は瞬間的な位相です。

   3. 瞬間エネルギー $$|a_i(t){|}^2$$、と瞬時周波数 $${\omega }_i(t)\equiv d{\theta }_i(t)/dt$$ を計算します。サンプルレートを指定した場合、hht は、$${\omega }_i(t)$$ を周波数 (Hz) に変換します。

   4. imfinse の瞬間エネルギーと imfinsf の瞬時周波数を出力します。

3. 出力引数なしで呼び出した場合、hht は、信号のエネルギーを、振幅に比例する色をもつ時間と周波数の関数としてプロットします。

参照

[1] Huang, Norden E, and Samuel S P Shen. Hilbert–Huang Transform and Its Applications. 2nd ed. Vol. 16. Interdisciplinary Mathematical Sciences. WORLD SCIENTIFIC, 2014. https://doi.org/10.1142/8804.

[2] Huang, Norden E., Zhaohua Wu, Steven R. Long, Kenneth C. Arnold, Xianyao Chen, and Karin Blank. “ON INSTANTANEOUS FREQUENCY.” Advances in Adaptive Data Analysis 01, no. 02 (April 2009): 177–229. https://doi.org/10.1142/S1793536909000096.

拡張機能

C/C++ コード生成

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

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

• 名前と値のペアを使用して指定する引数はコンパイル時の定数でなければなりません。

• コード生成では、timetable はサポートされていません。

参考

emd | hilbert | vmd

トピック

• 時間-周波数ギャラリー
• 解析信号とヒルベルト変換
• ヒルベルト変換および瞬時周波数

R2018a で導入

×

例を開く

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

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