hht

Hilbert-Huang transform

Syntax

hs = hht(imf,fs)
[hs,f,t] = hht(imf,fs)
[hs,f,t,imfinsf,imfinse] = hht(___)
[___] = hht(___,Name,Value)
hht(___)
hht(___,freqlocation)

Description

hs = hht(imf,fs) returns the Hilbert spectrum hs of the signal specified by intrinsic mode functions imf, and the optional sampling frequency fs. 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,f,t] = hht(imf,fs) returns frequency vector f and time vector t in addition to hs.

example

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

example

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

collapse all

For this example, consider 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.

Load the nonstationary signal data with the sampling frequency fs, and visualize the mixed sinusoidal signal.

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

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.

For this example, consider 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.

Load the nonstationary signal data with the sampling frequency fs, and visualize the mixed sinusoidal signal.

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

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 fs, 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.

Input Arguments

collapse all

Intrinsic mode function, specified as a matrix or timetable. 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. emd is used to decompose and simplify complicated signals into a finite number of intrinsic mode functions required to perform Hilbert spectral analysis.

hht treats each column in imf as an intrinsic mode function. For more information on computing imf, see emd.

Sampling frequency, specified as a positive scalar. Sampling frequency or sampling rate is the number of samples per unit time.

If fs is not supplied, a normalized frequency of is used to compute the Hilbert spectrum. If imf is specified as a timetable, the sampling time is inferred from it.

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 Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example: 'FrequencyResolution',1

Frequency limits to compute Hilbert spectrum, specified as the comma-separated pair consisting of 'FrequencyLimits' and a 1-by-2 integer-valued vector. FrequencyLimits is specified in Hz.

Frequency resolution to discretize frequency limits, specified as the comma-separated pair consisting of 'FrequencyResolution' and a positive scalar.

Specify FrequencyResolution 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, while f_low is the lower limit.

Minimum threshold value of Hilbert spectrum, specified as the comma-separated pair consisting of 'MinThreshold' and a scalar.

MinThreshold sets elements of hs to 0 when the corresponding elements of 10log10(hs) are less than MinThreshold.

Output Arguments

collapse all

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

Frequency vector of the signal, returned as a numeric array. 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:FrequencyResolution:f_high

Time vector of the signal, returned as a numeric array 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 imf is specified as an array.

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

Instantaneous frequency of each IMF, returned as an array or timetable.

insf has the same number of columns as imf and is returned as:

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

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

Instantaneous energy of each IMF, returned as an array or timetable.

inse has the same number of columns as imf and is returned as:

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

  • A timetable, if imf 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 decomposes the data set x into a finite number of intrinsic mode functions.

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

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

    2. Expresses zi as zi(t)=ai(t)ejθi(t), where ai(t) is the instantaneous amplitude and θi(t) is the instantaneous phase.

    3. Computes the instantaneous energy, |ai(t)|2, and the instantaneous frequency, ωi(t)dθi(t)/dt. If given a sample rate, hht converts ω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. Singapore: World Scientific, 2014.

[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. Vol. 1, No. 2, 2009, pp. 177–229.

Introduced in R2018a