hht

Hilbert-Huang transform

Syntax

hs = hht(IMFs)

hs = hht(IMFs,fs)

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

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

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

hht(___)

hht(___,freqlocation)

Description

hs = hht(IMFs) returns the Hilbert spectrum hs of the signal specified by intrinsic mode functions IMFs. hs is useful for analyzing signals that comprise a mixture of signals whose spectral content changes in time. Use hht to perform Hilbert spectral analysis on signals to identify localized features.

hs = hht(IMFs,fs) returns the Hilbert spectrum hs of a signal sampled at a rate fs.

[hs,f,t] = hht(___) returns frequency vector f and time vector t in addition to hs. These output arguments can be used with either of the previous input syntaxes.

[hs,f,t,imfinsf,imfinse] = hht(___) also returns the instantaneous frequencies imfinsf and the instantaneous energies imfinse of the intrinsic mode functions for signal diagnostics.

[___] = hht(___,Name=Value) estimates Hilbert spectrum parameters with additional options specified by one or more name-value arguments.

hht(___) with no output arguments plots the Hilbert spectrum in the current figure window. You can use this syntax with any of the input arguments in previous syntaxes.

hht(___,freqlocation) plots the Hilbert spectrum with the optional freqlocation argument to specify the location of the frequency axis. Frequency is represented on the y-axis by default.

Examples

Hilbert Spectrum of Quadratic Chirp

Generate a Gaussian-modulated quadratic chirp. Specify a sample rate of 2 kHz and a signal duration of 2 seconds.

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)

Use emd to visualize the intrinsic mode functions (IMFs) and the residual.

emd(q)

Compute the IMFs of the signal. Use the 'Display' name-value pair to output a table showing the number of sifting iterations, the relative tolerance, and the sifting stop criterion for each 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.

Use the computed IMFs to plot the Hilbert spectrum of the quadratic chirp. Restrict the frequency range from 0 Hz to 20 Hz.

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

Perform Empirical Mode Decomposition and Visualize Hilbert Spectrum of Signal

Load and visualize a nonstationary continuous signal composed of sinusoidal waves with a distinct change in frequency. The vibration of a jackhammer and the sound of fireworks are examples of nonstationary continuous signals. The signal is sampled at a rate fs.

load("sinusoidalSignalExampleData.mat","X","fs")
t = (0:length(X)-1)/fs;

plot(t,X)
xlabel("Time (s)")

The mixed signal contains sinusoidal waves with different amplitude and frequency values.

To create the Hilbert spectrum plot, you need the intrinsic mode functions (IMFs) of the signal. Perform empirical mode decomposition to compute the IMFs and residuals of the signal. Since the signal is not smooth, specify 'pchip' as the interpolation method.

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

The table generated in the command window indicates the number of sift iterations, the relative tolerance, and the sift stop criterion for each generated IMF. This information is also contained in info. You can hide the table by adding the 'Display',0 name value pair.

Create the Hilbert spectrum plot using the imf components obtained using empirical mode decomposition.

hht(imf,fs)

The frequency versus time plot is a sparse plot with a vertical color bar indicating the instantaneous energy at each point in the IMF. The plot represents the instantaneous frequency spectrum of each component decomposed from the original mixed signal. Three IMFs appear in the plot with a distinct change in frequency at 1 second.

Hilbert Spectrum of Whale Song

Load a file that contains audio data from a Pacific blue whale, sampled at 4 kHz. The file is from the library of animal vocalizations maintained by the Cornell University Bioacoustics Research Program. The time scale in the data is compressed by a factor of 10 to raise the pitch and make the calls more audible. Convert the signal to a MATLAB® timetable and plot it. Four features stand out from the noise in the signal. The first is known as a trill, and the other three are known as moans.

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

Use emd to visualize the first three intrinsic mode functions (IMFs) and the residual.

emd(whale,'MaxNumIMF',3)

Compute the first three IMFs of the signal. Use the 'Display' name-value pair to output a table showing the number of sifting iterations, the relative tolerance, and the sifting stop criterion for each 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.

Use the computed IMFs to plot the Hilbert spectrum of the signal. Restrict the frequency range from 0 Hz to 1400 Hz.

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

Compute the Hilbert spectrum for the same range of frequencies. Visualize the Hilbert spectra of the trill and moans as a mesh plot.

[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')

Compute Hilbert Spectrum Parameters of Signal

Load and visualize a nonstationary continuous signal composed of sinusoidal waves with a distinct change in frequency. The vibration of a jackhammer and the sound of fireworks are examples of nonstationary continuous signals. The signal is sampled at a rate fs.

load("sinusoidalSignalExampleData.mat","X","fs")
t = (0:length(X)-1)/fs;

plot(t,X)
xlabel("Time (s)")

The mixed signal contains sinusoidal waves with different amplitude and frequency values.

To compute the Hilbert spectrum parameters, you need the IMFs of the signal. Perform empirical mode decomposition to compute the intrinsic mode functions and residuals of the signal. Since the signal is not smooth, specify 'pchip' as the interpolation method.

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

The table generated in the command window indicates the number of sift iterations, the relative tolerance, and the sift stop criterion for each generated IMF. This information is also contained in info. You can hide the table by specifying 'Display' as 0.

Compute the Hilbert spectrum parameters: Hilbert spectrum hs, frequency vector f, time vector t, instantaneous frequency imfinsf, and instantaneous energy imfinse.

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

Use the computed Hilbert spectrum parameters for time-frequency analysis and signal diagnostics.

VMD of Multicomponent Signal

Generate a multicomponent signal consisting of three sinusoids of frequencies 2 Hz, 10 Hz, and 30 Hz. The sinusoids are sampled at 1 kHz for 2 seconds. Embed the signal in white Gaussian noise of variance 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));

Compute the IMFs of the noisy signal and visualize them in a 3-D plot.

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')

Use the computed IMFs to plot the Hilbert spectrum of the multicomponent signal. Restrict the frequency range to [0, 40] Hz.

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

Compute Hilbert Spectrum of Vibration Signal

Simulate a vibration signal from a damaged bearing. Compute the Hilbert spectrum of this signal and look for defects.

A bearing with a pitch diameter of 12 cm has eight rolling elements. Each rolling element has a diameter of 2 cm. The outer race remains stationary as the inner race is driven at 25 cycles per second. An accelerometer samples the bearing vibrations at 10 kHz.

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

The vibration signal from the healthy bearing includes several orders of the driving frequency.

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;

A resonance is excited in the bearing vibration halfway through the measurement process.

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

The resonance introduces a defect in the outer race of the bearing that results in progressive wear. The defect causes a series of impacts that recur at the ball pass frequency outer race (BPFO) of the bearing:

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

where $$f_0$$ is the driving rate, $$n$$ is the number of rolling elements, $$d$$ is the diameter of the rolling elements, $$p$$ is the pitch diameter of the bearing, and $$\theta$$ is the bearing contact angle. Assume a contact angle of 15° and compute the BPFO.

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

Use the pulstran function to model the impacts as a periodic train of 5-millisecond sinusoids. Each 3 kHz sinusoid is windowed by a flat top window. Use a power law to introduce progressive wear in the bearing vibration signal.

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 = 49e3;
nSamples = 1e5;
yImpact = pulstran(t,tx',xImpact,fs)/5;
yImpact = [zeros(1,nWear) yImpact(1,(nWear+1):nSamples)];

Generate the BPFO vibration signal by adding the impacts to the healthy bearing signal. Plot the signal and select a 0.3-second interval starting at 5.0 seconds.

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

Zoom in on the selected interval to visualize the effect of the impacts.

xlim([xLimLeft xLimRight])

Add white Gaussian noise to the signals. Specify a noise variance of $$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")

Use emd to perform an empirical mode decomposition of the healthy bearing signal. Compute the first five intrinsic mode functions (IMFs). Use the 'Display' name-value argument to output a table showing the number of sifting iterations, the relative tolerance, and the sifting stop criterion for each 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.

Use emd without output arguments to visualize the first three IMFs and the residual.

emd(yGood,MaxNumIMF=5)

Compute and visualize the IMFs of the defective bearing signal. The first empirical mode reveals the high-frequency impacts. This high-frequency mode increases in energy as the wear progresses.

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)

Plot the Hilbert spectrum of the first empirical mode of the defective bearing signal. The first mode captures the effect of high-frequency impacts. The energy of the impacts increases as the bearing wear progresses.

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

The Hilbert spectrum of the third mode shows the resonance in the vibration signal. Restrict the frequency range from 0 Hz to 100 Hz.

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

For comparison, plot the Hilbert spectra of the first and third modes of the healthy bearing signal.

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

Input Arguments

IMFs — Intrinsic mode functions

matrix | timetable

Intrinsic mode functions, specified as a matrix or timetable. The hht function treats each column in IMFs as an intrinsic mode function (IMF). An IMF is any signal whose envelope is symmetric with respect to zero and whose numbers of extrema and zero crossings differ by at most one.

Use emd (empirical mode decomposition) or vmd (variational mode decomposition) to decompose and simplify complicated signals into a finite number of IMFs required to perform Hilbert spectral analysis.

fs — Sample Rate

2π (default) | positive scalar

Sample rate, specified as a positive scalar. If fs is not supplied, a normalized frequency of 2π is used to compute the Hilbert spectrum. If IMFs is specified as a timetable, the sample rate is inferred from it.

freqlocation — Location of frequency axis on plot

"yaxis" (default) | "xaxis"

Location of frequency axis on the plot, specified as "yaxis" or "xaxis". To display frequency data on the y-axis or x-axis of the plot, specify freqlocation as "yaxis" or "xaxis" respectively.

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: 'FrequencyResolution',1

FrequencyLimits — Frequency limits to compute Hilbert spectrum

[0,fs/2] (default) | 1-by-2 integer-valued vector

Frequency limits to compute Hilbert spectrum, specified as a 1-by-2 integer-valued vector. FrequencyLimits is specified in Hz.

FrequencyResolution — Frequency resolution to discretize frequency range

(f_high-f_low)/100 (default) | positive scalar

Frequency resolution to discretize frequency limits, specified as a positive scalar. FrequencyResolution is specified in Hz. If FrequencyResolution is not specified, a value of (f_high-f_low)/100 is inferred from FrequencyLimits. Here, f_high is the upper limit of FrequencyLimits and f_low is the lower limit.

MinThreshold — Minimum threshold value of Hilbert spectrum

-inf (default) | scalar

Minimum threshold value of Hilbert spectrum, specified as a scalar. MinThreshold sets elements of hs to 0 when the corresponding elements of $$10{\log }_{10}\left(hs\right)$$ are less than MinThreshold.

Output Arguments

hs — Hilbert spectrum of signal

sparse matrix

Hilbert spectrum of the signal, returned as a sparse matrix. Use hs for time-frequency analysis and to identify localized features in the signal.

f — Frequency values

vector

Frequency values of the signal, returned as a vector. hht uses the frequency vector f and the time vector t to create the Hilbert spectrum plot.

Mathematically, f is denoted as: f = f_low : f_res : f_high, where f_res is the frequency resolution.

t — Time values

vector | duration array

Time values of the signal, returned as a vector or a duration array. hht uses the time vector t and the frequency vector f to create the Hilbert spectrum plot.

t is returned as:

• An array, if IMFs is specified as an array.

• A duration array, if IMFs is specified as a uniformly sampled timetable.

imfinsf — Instantaneous frequency of each IMF

vector | matrix | timetable

Instantaneous frequency of each IMF, returned as a vector, a matrix, or a timetable.

imfinsf has the same number of columns as IMFs and is returned as:

• A vector, if IMFs is specified as a vector.

• A matrix, if IMFs is specified as a matrix.

• A timetable, if IMFs is specified as a uniformly sampled timetable.

imfinse — Instantaneous energy of each IMF

vector | matrix | timetable

Instantaneous energy of each IMF, returned as a vector, a matrix, or a timetable.

imfinse has the same number of columns as IMFs and is returned as:

• A vector, if IMFs is specified as a vector.

• A matrix, if IMFs is specified as a matrix.

• A timetable, if IMFs is specified as a uniformly sampled timetable.

Algorithms

The Hilbert-Huang transform is useful for performing time-frequency analysis of nonstationary and nonlinear data. The Hilbert-Huang procedure consists of the following steps:

1. emd or vmd decomposes the data set x into a finite number of intrinsic mode functions.

2. For each intrinsic mode function, x_i, the function hht:

   1. Uses hilbert to compute the analytic signal, $$z_i(t)=x_i(t)+jH\{x_i(t)\}$$, where H{x_i} is the Hilbert transform of x_i.

   2. Expresses z_i as $$z_i(t)=a_i(t)e^{j{\theta }_i(t)}$$, where a_i(t) is the instantaneous amplitude and $${\theta }_i(t)$$ is the instantaneous phase.

   3. Computes the instantaneous energy, $$|a_i(t){|}^2$$, and the instantaneous frequency, $${\omega }_i(t)\equiv d{\theta }_i(t)/dt$$. If given a sample rate, hht converts $${\omega }_i(t)$$ to a frequency in Hz.

   4. Outputs the instantaneous energy in imfinse and the instantaneous frequency in imfinsf.

3. When called with no output arguments, hht plots the energy of the signal as a function of time and frequency, with color proportional to amplitude.

References

[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.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• Name-value arguments must be compile-time constants.

• Due to the nature of the binning algorithm used in the estimation of the instantaneous frequencies, there may be a mismatch between results obtained using MATLAB^® code and generated C/C++ code. The mismatches are due to off-by-one differences in the frequency estimation.

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

The hht function supports GPU array input with these usage notes and limitations:

• Due to the nature of the binning algorithm used in the estimation of the instantaneous frequencies, there may be a mismatch between results obtained on a CPU and GPU. The mismatches are due to off-by-one differences in the frequency estimation.

For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced in R2018a

R2024b: Use gpuArray objects

The hht function supports gpuArray objects. You must have Parallel Computing Toolbox™ to use this functionality.