このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
hht
ヒルベルト・ファン変換
構文
説明
[___] = hht(___,
では、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) で繰り返される一連の影響を引き起こします。
ここで、 は駆動レート、 は回転要素の数、 は回転要素の直径、 はベアリングのピッチの直径、 はベアリングの接触角です。接触角は 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])
ホワイト ガウス ノイズを信号に付加します。 のノイズ分散を指定します。
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.017132 | SiftMaxRelativeTolerance 2 | 3 | 0.12694 | SiftMaxRelativeTolerance 3 | 6 | 0.14582 | SiftMaxRelativeTolerance 4 | 1 | 0.011082 | SiftMaxRelativeTolerance 5 | 2 | 0.03463 | SiftMaxRelativeTolerance Decomposition stopped because maximum number of intrinsic mode functions was extracted.
emd
を出力引数なしで使用し、最初の 3 つの IMF と残差を可視化します。
emd(yGood,MaxNumIMF=5)
欠陥のあるベアリング信号の IMF を計算および可視化します。最初の経験的モードでは、高周波数に影響が見られます。この高周波数モードでは、摩耗が進行するにつれてエネルギーが増加します。
imfBad = emd(yBad,MaxNumIMF=5,Display=1);
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 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])
入力引数
fs
— サンプル レート
2π
(既定値) | 正のスカラー
サンプル レート。正のスカラーで指定します。fs
が指定されない場合、2π
の正規化された周波数がヒルベルト スペクトルの計算に使用されます。imf
が timetable として指定されている場合、これからサンプル レートが推定されます。
freqlocation
— プロット上の周波数軸の位置
"yaxis"
(既定値) | "xaxis"
プロット上の周波数軸の位置。"yaxis"
または "xaxis"
として指定します。プロットの y 軸または x 軸上に周波数データを表示するには、それぞれ "yaxis"
または "xaxis"
として freqlocation
を指定します。
名前と値の引数
オプションの引数のペアを Name1=Value1,...,NameN=ValueN
として指定します。ここで、Name
は引数名で、Value
は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。
R2021a より前では、コンマを使用して名前と値をそれぞれ区切り、Name
を引用符で囲みます。
例: 'FrequencyResolution',1
FrequencyLimits
— ヒルベルト スペクトルを計算する周波数限界
[0,fs
/2]
(既定値) | 1 行 2 列の整数値のベクトル
fs
/2]ヒルベルト スペクトルを計算する周波数限界。1 行 2 列の整数値ベクトルとして指定します。FrequencyLimits
は Hz で指定します。
FrequencyResolution
— 周波数範囲を離散化する周波数分解能
(f_high-f_low)/100 (既定値) | 正のスカラー
周波数限界を離散化する周波数分解能。正のスカラーとして指定します。FrequencyResolution
は Hz で指定します。FrequencyResolution
が指定されない場合、(fhigh-flow)/100 の値は FrequencyLimits
から推定されます。ここで、fhigh は FrequencyLimits
の上限、flow は下限です。
MinThreshold
— ヒルベルト スペクトルの最小しきい値
-inf
(既定値) | スカラー
ヒルベルト スペクトルの最小しきい値。スカラーとして指定します。MinThreshold
は、 の対応する要素が MinThreshold
より小さい場合に hs
の要素を 0 に設定します。
出力引数
hs
— 信号のヒルベルト スペクトル
スパース行列
信号のヒルベルト スペクトル。スパース行列として返されます。hs
を使用して時間-周波数を解析し、信号の局在的な特徴を特定します。
f
— 周波数値
ベクトル
信号の周波数値。ベクトルとして返されます。hht
は、周波数ベクトル f
と時間ベクトル t
を使用してヒルベルト スペクトル プロットを作成します。
数学的には、f
は f = flow : fres : fhigh と表されます。ここで、fres は周波数分解能です。
imfinsf
— 各 IMF の瞬時周波数
ベクトル | 行列 | timetable
各 IMF の瞬時周波数。ベクトル、行列、または timetable として返されます。
imfinsf
は imf
と同じ列数であり、以下として返されます。
imf
をベクトルとして指定した場合はベクトル。imf
を行列として指定した場合は行列。imf
を等間隔サンプルの timetable として指定した場合は timetable。
imfinse
— 各 IMF の瞬間エネルギー
ベクトル | 行列 | timetable
各 IMF の瞬間エネルギー。ベクトル、行列、または timetable として返されます。
imfinse
は imf
と同じ列数であり、以下として返されます。
imf
をベクトルとして指定した場合はベクトル。imf
を行列として指定した場合は行列。imf
を等間隔サンプルの timetable として指定した場合は timetable。
アルゴリズム
ヒルベルト・ファン変換は、非定常および非線形データの時間-周波数解析の実行に役立ちます。ヒルベルト・ファン手順は以下のステップから構成されます。
各固有モード関数 xi の場合、関数
hht
は次のとおりです。出力引数なしで呼び出した場合、
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++ コードを生成します。
使用上の注意および制限:
名前と値の引数はコンパイル時の定数でなければなりません。
バージョン履歴
R2018a で導入
MATLAB コマンド
次の MATLAB コマンドに対応するリンクがクリックされました。
コマンドを MATLAB コマンド ウィンドウに入力して実行してください。Web ブラウザーは MATLAB コマンドをサポートしていません。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)