Hilbert-Huang Transform

This submission is a realization of the Hilbert-Huang transform (HHT).

Description

The function plot_hht is a realization of the Hilbert-Huang transform (HHT). The HHT decomposes a signal into intrinsic mode functions (or IMFs), and obtain the instantaneous frequency data. It is designed to work well for data that are nonstationary and nonlinear (http://en.wikipedia.org/wiki/Hilbert-Huang_Transform). Learn more about the HHT from the attached pdf.

The function plot_hht is best used with spsnip_gui available at http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=11002&objectType=FILE.

Essential files: plot_hht.m, emd.m, findpeaks.m
Accompanying files: Hum.wav, HHT.pdf

Acknowledgements

This file inspired Imf For Bearing Fault Diagnosis.

• Required Products: Signal Processing Toolbox; Spline Toolbox
• MATLAB release: MATLAB 5.3.1 (R11.1)

[junqua]

 

[RAM KINKER MISHRA]

Alan Tan, I had a query in emd decomposition you have assigned 'SD=inf' but in my opinion instead of that sd should be calculated by using intermediate IMF from starting onwards.

[Jaidev Deshpande]

Hi Alan,

I have proposed to come up with a Python implementation of the HHT under the Google Summer of Code programme. But I'm worried about the patented status of the HHT.

I see many implementations of the HHT that are freely available, and a lot of research papers are based on HHT too. If the interest is purely academic and non-commercial, can HHT be used freely?

Please help me out.

Thanks

[dedenlaut]

How to plot this result in energy-frequency-time domain?

[Sergi Andruschenko]

 

[RAM KINKER MISHRA]

This is the code i am using.

function hht_meanfreq(x,Ts)
imf = emd(xsin);
%emd finds the decomposed segments
for i=1:length(imf)
imf1=cell2mat(imf(i));
%conversion of cell to matrix as i am %use to work with matrix
a1(i,:)=imf1;
a(i,:)=abs(hilbert(a1(i,:)));
%absolute value of each hilbert %transformed IMF
th(i,:)=angle(hilbert(a1(i,:)));
%fINDING PHASE OF EACH IMF
d(i,:) = diff(th(i,:))/Ts/(2*pi);
%CALCULAION OF INSTANTANEOUS FREQUENCY
end

for j=1:length(imf)
mnf1(j)
=sum(d(j,:).*(a(j,1:end-1)).^2)/norm(a(j,1:end-1))^2;
%CALCULATION OF MEAN %INSTATNTANEOUS %FREQUENCY OF EACH IMF
end
for z=1:length(imf)
mnf=sum(norm(a(z,:).*mnf1(z)))/sum(norm(a(z,:)));
%FINDING MEAN FREQUENCY
end

[MANOJ KUMAR PALO]

please provide me the matlab code of hilbert -transform

[RAM KINKER MISHRA]

My aim is to find the mean frequency derived via Hilbert-Huang transform.
Algorithm
1. first find the mean instantaneous frequency for each IMF.
2.Then using these instantaneous frequencies i will find out the mean frequency of the overall signal.
I have tried to implement the algorithm in the above program but i am not able to get the result properly.

[RAM KINKER MISHRA]

function hht_meanfreq(x,Ts)
imf = emd(xsin);
%emd finds the decomposed segments
for i=1:length(imf)
imf1=cell2mat(imf(i));
%conversion of cell to matrix as i am use to work with matrix
a1(i,:)=imf1;
a(i,:)=abs(hilbert(a1(i,:)));
% absolute value of each hilbert transformed IMF
th(i,:)=angle(hilbert(a1(i,:)));
%fINDING PHASE OF EACH IMF
d(i,:) = diff(th(i,:))/Ts/(2*pi);
%CALCULAION OF INSTANTANEOUS FREQUENCY
end

for j=1:length(imf)
mnf1(j)=sum(d(j,:).*(a(j,1:end-1)).^2)/norm(a(j,1:end-1))^2;
%CALCULATION OF MEAN INSTATNTANEOUS FREQUENCY OF EACH IMF
end
for z=1:length(imf)
mnf=sum(norm(a(z,:).*mnf1(z)))/sum(norm(a(z,:)));
%FINDING MEAN FREQUENCY
end

[Ali Obodi]

To Author,

when i run plot_hht.m, i had following message appeared on command window:

??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N) to change the limit. Be aware that exceeding your available stack space cancrash MATLAB and/or your computer.

Error in ==> fileparts

Then, i used
set(0,'RecursionLimit',1000) (maximum value 1000) and computer crashes every time! when i used lower values i had the same previous message.
Please need advice regarding this problem. Many thanks

[pankaj badoni]

hi please anyone tell how to apply this one on images esp of medical images.

plz..reply fast.

[aami]

cud any one plz give me the link to the code for decomposing an image using EMD...
Thanx in advance...

Meenu

[Lukasz]

Shouldn't you unwrap the angle th before sending it to diff?

[Omar]

I'm interested in using the envelope of my signal that might I get using your code for further analysis. But, by plotting both original signal and its envelope on the same figure, it does not show me that the envelope follow the peaks of my original signal as shown in figure 1.2 (upper envelope) in the HHT.pdf file. So, could you please suggest me an idea how to make the envelope follow the local maxima.

Thanks again,
Omar

[Omar]

 

[Bubo Bubic]

Hi,
I got next message trying to run plot_hht.m:

??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N) to change the limit. Be aware that exceeding your available stack space cancrash MATLAB and/or your computer.

Error in ==> fileparts

Then I run from command line
set(0,'RecursionLimit',2600) (maximum value 2600) and computer crashes every time! For lower values got the same above message.
Please help with suggestion or advices.
Thanks

[Paul Harfas]

Hi.
First, I must say this piece of code helped me with my work, thank you for your effort.
I am currently working on improving your code by adding extrema at the beginning and end of the data series based on the actual data points (slope-based method), not just adding zeros. Adding zeros will result in "fake" (artificially induced) oscillations that will propagate through adjacent modes (IMFs) and induce errors. I will post back when I reach a result.

@rakesh: to obtain a surface (3D) plot of the time-frequency-energy distribution, you should plot abs(hilbert(imf{k})).^2 as the value on the Z axis (plot_hht.m line 15).

[Alan Tan]

@Rakesh: Not sure what you meant with the 3d plot and vertical axis. Sorry.
@Hans: I have not experimented with files of .edf extension; I can't help much with that -- sorry.

[Hans]

To Author,

I have another one more question is that, my data is .edf but once i plot the hht using hht_plot.m, the result came out is empty but IMF is running and not empty, is there something wrong with my data?

[Rakesh]

The 1st two plots are the plots of 2 imfs with highest engery content. However, according to the theory, all imfs (and hence their frqs) are plotted with respect to time. I did "hold on" and have plotted time-freq distribution for all the imfs. the plot looks kind of crowded but still gives enough information. However i still hav following questions,

1. in hilbert spectrum, in 3D plot, vertical axis amplitude is function of time. but u have obtained a single value of energy for each imf. but actually it should have been function of time.Am i right?

2. How should I get Spactrogram type of plot?

[Mark Shore]

Just tried this, no comments yet except that the spline toolbox is not required (at least for 2009b). A call is made to spline.m which is found in MATLAB's standard toolbox in the polyfun directory.

[Alan Tan]

Hi Jaime, thank you for the comment.

From Matlab's help on the supplied hilbert.m function, it appears that there is no normalization involved. I had no particular reason on the choice of the figure colours -- you can always change the colour of the plot by modifying Line 25 of the plot_hht.m function.

[Jaime Delgado Saa]

Hello, Thanks for this tool, It works fine for me.One question: do you use normalized hilbert transform for the calculation of the instantaneous frequency of the imf's.

just curious : why the pictures are black and white?

[Alan Tan]

The following reply is given to Agnes over email.

Hi Agnes, The 2 initial plots are the time-frequency plots of the HHT corresponding to the 2 IMFs with of the highest energy. Subsequent plots are the IMFs, plotted against time.

[Agnes]

Please describe me what the plots are exactelly. I suppose that the first is the HHT of the datas, but what is the second?what is the difference between them? And the other plots are the various envelopes?
And if I have a 0.7 KHz sampling rate, than I have to give 700 as the Ts or 1/700?
Thanks for your help!

[Raymond Cheng]

Thanks for your sharing.

[delinda he]

Would you guys please help me to figure out how does the Hilbert transform work??
function plot_hht(x,Ts)
% Plot the HHT.
% plot_hht(x,Ts)
% :: Syntax
% The array x is the input signal and Ts is the sampling period.
% Example on use: [x,Fs] = wavread('Hum.wav');
% plot_hht(x(1:6000),1/Fs);
% Func : emd

% Get HHT.
imf = emd(x); %emd is a function, it returns the IMF (intrinsic mode functions)

%%%%%%% The following one is to calculate Hilbert transform for each IMFs, but I do no know how it works???????@#$%%

for k = 1:length(imf)
b(k) = sum(imf{k}.*imf{k});
th = angle(hilbert(imf{k}));
d{k} = diff(th)/Ts/(2*pi);
end
[u,v] = sort(-b);
b = 1-b/max(b);

% Set time-frequency plots. I do not how it plot the time-frequency figures too
N = length(x);
c = linspace(0,(N-2)*Ts,N-1);
for k = v(1:2)
figure, plot(c,d{k},'k.','Color',b([k k k]),'MarkerSize',3);
set(gca,'FontSize',8,'XLim',[0 c(end)],'YLim',[0 1/2/Ts]); xlabel('Time'), ylabel('Frequency');
end

Thanks for your help.

[Alan Tan]

For the benefit of other readers, the following is the reply given to Yu over email.

Hi Yu, I am assuming the s you mentioned is output of the getspline function. The array s is formed by invoking Matlab's spline function, i.e., spline([0 p N+1],[0 x(p) 0],1:N). You'd note that additional endpoint conditions (loosely put, the spline should settle towards 0 at the endpoints) have been imposed as a means to contain the endpoints (of the resulting spline) to within reasonable bounds. Removing these conditions may result in the 'extreme peaks' you mentioned. I hope that clears your doubt.

[Yu]

to Auther,

hi, you did a very good job, and with your function i saved a lot of time. i have a question about thd 's'- sensitivity. i got different extreme peaks , when i change the 's'. can you tell me whats the meaning of 's'?
thank you

[Alan Tan]

Hi Asim, Ts is a single number (not a column of numbers) denoting the sampling time. If you have samples collected every, say 1 day, then set Ts = 1. Please send me an email if you have other doubts -- do not post your questions here since I do not monitor this site regularly. Thank you.

[Asim]

To Author, can you please explain how we are defining the variable 'Ts'? If I have a dataseries collected everyday for 45 years, will this column 1 will be Ts or just the value 1? In my dataset I have two column with second with dataset and first with day of measurement (regular basis).

[SQ Yuan]

 

[Kechang Fu]

 

[Legolas Long]

why can't i use it correctly?

[Chian Wong]

The user would also need the Signal Processing Toolbox to perform the Hilbert Transform itself, though this is elementary to code.

[The Author]

If you don't have phase.m in your version of Matlab (I have in mine!), you can substitute with angle.m in the plot_hht function.

[The Author]

Forgot to add that you will need to have the Spline Toolbox to run this one.