Hilbert-Huang Transform: matlab Hilbert-Huang Transform: matlabの実装

Hilbert-Huangのmatlab実装、資料のまとめ、割と雑多です...。Web上の投稿者の方々に感謝します :)

コア : 次のコードは、HHT のマージナルスペクトルとそれに対応する周波数を計算します。
ツールキットの要件 : G-Rilling EMD Toolbox、TFTB Toolbox
アタッチメント 黄Eさんのグループが開発したツールキット（下記からダウンロードできます。 こちら が見つかりました）、ここでは使用しません。

% Empirical mode decomposition, resulting in intrinc mode functions.
% Without parameter 'MAXMODES' the process may be seriously delayed by
% decompose original signals into too many IMFs (not necessary, 9 is 
% enough generally)
imfs = emd(oriSig, 'MAXMODES', 9);
% HHT spectrum: hhtS
[A, f, t] = hhspectrum(imfs);
[hhtS, ~, fCent] = toimage(A, f, t);
% Marginal hilbert spectrum: hhtMS, xf: correspondig frequency
for k = 1 : size(hhtS, 1)
    hhtMS(k) = sum(hhtS(k, : )) * 1 / fs;
end
xf = fCent(1, :) . * fs;

G-Rilling EMD Toolbox Toolkit関連設定: 咗凡春秋

簡単に言うと、パスを設定した後、入力します。 インストール で完了です。

EMD分解後のヒルベルトスペクトルのMatlab解析（相関関数はすべてG-Rilling EMD Toolboxに帰属します。）

• hhspectrumの機能説明 (8階：旧旧学生)

% [A,f,tt] = HHSPECTRUM(imf,t,l,aff)
% Input:
%- imf : matrix with one IMF per row % Substitute the IMF obtained from the emd decomposition, which is the c-variable written in your program, without adding the trend term in the last row
%- t : time instants % instantaneous time or duration ? (write [1: signal length] on it, the real time can be converted according to the sampling rate)
%- l : estimation parameter for instfreq % The estimation parameter for the instantaneous frequency ? (write 1 to determine the boundary from the first point in the estimation of the instantaneous frequency)
%- aff : if 1, displays the computation evolution % option to display the computation process, write 0 if you don't want to display it
%
% Output:
% - A : amplitudes of IMF rows
% - f : instantaneous frequencies
% - tt : truncated time instants % Cutoff time ? (the truncated time, which returns the time corresponding to the instantaneous frequencies, should be shorter than the original signal time, determined by the previous l value)
% calls
% calls:
% - hilbert : computes the analytic signal
% - instfreq : computes the instantaneous frequency % instantaneous frequency
%
% Example.
[A,f,tt] = hhspectrum(c(1:end-1,:),[1:n],1,0);
[im,tt,ff] = toimage(A,f,tt,512);
disp_hhs(im,tt,[],fs);

• 9階：旧旧学生

時間周波数図の概念については、ウェーブレットやガボールなどの時間周波数解析の共同手法と連動していると思います。

ウェーブレットやガボールなどはマルチスケール解析の概念を持ち、得られる時間-周波数分布は2次元行列（横軸が時間、縦軸が周波数で、異なる色（スペクトログラム）や大きさの異なるウォーターフォールプロットで表現できる）です。

HHTの場合、気体-時間周波数プロットの概念がやや異なると思います。EMD分解の役割は、瞬時周波数の概念を適用できるように複雑な信号を単純な1成分信号に離散化することであり、ヒルベルト変換の目的は瞬時周波数を分析することです。つまり、HHTはウェーブレットなどのような一連の値（マルチスケール解析）ではなく、各瞬間に1つの値だけを取得するわけです。ですから、その時間-周波数分布から見えるのは、グラフではなく線なのです。

• 27階：ソングジー41
  HHTの重要な点は、EMDで分解した後、IMPをヒルベルトとして、各固有モード関数だけが過渡的な情報を持つようにすることである。原信号は多くの固有モード関数の和であり、全ての原信号に対して過渡的な情報は存在しない。これがHHTの基本的な出発点です。

G-Rilling EMD Toolboxに基づくヒルベルトスペクトル計算のサンプルコード: nike815

% Example program
N=1000;
n=1:N;
fs=1000;
t=n/fs;
fx=50;
fy=150;
x=2*sin(2*pi*fx*t);
y=sin(2*pi*fy*t);
z=x+y;
data=z;
imf=emd(data); % perform EMD decomposition on the input signal    
[A,f,t]=hhspectrum(imf); % instantaneous frequency and amplitude for IMF components: A: is the amplitude vector of each IMF, f: the instantaneous frequency corresponding to each IMF, t: time series number
[E,t,Cenf]=toimage(A,f); % synthesize each IMF signal to obtain the Hilbert spectrum, E: the corresponding amplitude value, Cenf: the corresponding center frequency of each grid. Here the horizontal axis is time and the vertical axis is frequency. Immediate frequency plot (with color indicating the magnitude of the third dimensional value) and 3D plot (3D coordinate system: time, central frequency, amplitude)        
cemd_visu(data,1:length(data),imf); % display each IMF component and residual signal
disp_hhs(E); % Hilbert spectrum
% draw the marginal spectrum
colormap(flipud(gray)); % display in black and white
%N=length(Cenf); % Set the number of frequency points. Completely from the theoretical formula. The center frequency is important after gridding, and we all think about it from the perspective of continuous data becoming discrete, I believe it should be easy to understand
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;
end
figure(3);
plot(Cenf(1,:)*fs,bjp); % Make a marginal spectrum plot The frequencies have been sampled into discrete frequencies with a certain window length when the Hilbert spectrum is obtained, so the frequency axis is already the center frequency at this point
xlabel('frequency / Hz');
ylabel('amplitude');

% Plot the instantaneous envelope and instantaneous frequency graph
figure;
subplot(221),plot(t/N,A(1,:));xlabel('time t/s');ylabel('amplitude');title('imf1 component transient envelope');
subplot(222),plot(t/N,f(1,:)*fs);xlabel('time t/s');ylabel('frequency');title('imf1 component transient frequency');
subplot(223),plot(t/N,A(2,:));xlabel('time t/s');ylabel('amplitude');title('imf2 component transient envelope');
subplot(224),plot(t/N,f(2,:)*fs);xlabel('time t/s');ylabel('frequency');title('imf2 component transient frequency');

ヒルベルト変換-ヒルベルトスペクトル、ヒルベルト周辺スペクトル、包絡線スペクトル。兪慶民

ヒルベルトスペクトル : fftに続く信号のヒルベルト変換で、全周波数帯域における信号振幅の時間および周波数によるパターンを表す。
ヒルベルト限界スペクトル ヒルベルト限界スペクトル：ヒルベルトスペクトルの積分で、全周波数帯域における信号振幅の周波数に対する変化を表し、フーリエスペクトルと同等であるが、フーリエスペクトルよりも周波数分解能が高いのが特徴。
ヒルベルト包絡線スペクトル ヒルベルト変換。Hilbert変換した包絡線をFFTし、HilbertスペクトルやHilbert限界スペクトルとは異なり、直接Hilbert変換した信号で解析関数を構築し、解析関数に基づいてモード値を求め、モード値は包絡線であり、信号包絡線をFFTしてHilbert包絡線スペクトルを取得することである。

学習メモ：ローゲン

限界スペクトルは、統計的な意味で、データ全体の各周波数点における累積振幅分布を表しているが、フーリエスペクトルのある点における振幅は、信号全体にその周波数を含むデルタ関数成分が存在することを示し、振幅が大きいと、単にデータセグメント全体にわたって局所的に存在する可能性が高いことを表している。
正確な時間を知りたいのであれば、時間-振幅-周波数の3次元スペクトルであるHHTスペクトルを見ればいいのです。瞬間的な周波数に関して言えば、フーリエ変換は局在化を強調するのではなく、グローバル化を強調する。瞬間周波数の独自の定義を提案しているのは、私たちのHHTなのです。したがって、フーリエ変換の測定にも瞬時周波数をとるのはフェアではない。
限界スペクトルについては、ヒルベルトスペクトルの時間的な積分であり、積分の観点からは、任意の1つの周波数の次数について時間的にすべての大きさを加算することと同じであり、この周波数の次数についてすべての時間の大きさの累積を反映したものである 周波数の振幅は、その周波数が全時間に出現する可能性と表現されていますが、個人的には、出現するのだから存在する、可能性とは言えない、ただ、出現した回数を累積したものが振幅になると考えており、フーリエの場合は、そのように考えています。フーリエの場合は、サイン・コサイン関数で信号をフィッティングするときに、この次数の周波数が必要だと説明できるだけで、振幅が大きいということは、信号のフィッティングに大きく貢献しているということで、発生の可能性を示すものではありません。
エネルギースペクトルは、限界スペクトルの2乗と理解できる。これは単にエネルギーの形をしているだけで、正確にエネルギーを表しているかどうかは、さらに調べる必要がある。また、瞬時エネルギーがあり、これらの用語は本当に魅力的だが、具体的にどのように、解決する必要がある。

より詳細なHHTのソース（matlabコード付き）: cwjy

• マージナルスペクトルとフーリエスペクトルの比較
  限界スペクトルはデータ全体の各周波数点での累積振幅分布を統計的に特徴づけるものであり、フーリエスペクトルのある点周波数での振幅は、信号全体にその周波数を含むデルタ関数成分が存在することを示すもので、意味が異なる。
  異なる効果 限界スペクトルは、非平滑な信号を扱うことができます。信号にある周波数のエネルギーが現れているということは、その周波数の振動波が現れているはずで、つまり限界スペクトルは、信号の実際の周波数成分をより正確に反映させることができるのです。そしてフーリエ変換は、滑らかな信号しか扱えない。
  
• HHT & ヒルベルト変換
  ヒルベルト変換は、単に信号の瞬間的な振幅、周波数、位相を求めるもので、意味のない負の周波数が存在する可能性があります。HHT変換では、まず信号をEMDに分解し、異なるスケールでの成分を与え、各成分のヒルベルト変換により、本当の意味を持つ瞬間的な周波数を求めます。