This submission is a realization of the HilbertHuang transform (HHT).
The function plot_hht is a realization of the HilbertHuang 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/HilbertHuang_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
1.1  Added BSD license. 

Replaced phase() with angle() in plot_hht(). 

Added requirement for the Spline Toolbox. 
Inspired: IMF for Bearing Fault Diagnosis
Mingjun Xiang (view profile)
This comment system is stupid!I can rate infinate times!
Mingjun Xiang (view profile)
Now,my matlab's version is not well worked with this code.I need to change 'wavread' to 'audioread',and I find that matlab has define 'findpeaks',so if i move this code to one .m,It is some wrong.
Mingjun Xiang (view profile)
Well,It gaves me result!Thank u!
James Too (view profile)
Li Liang (view profile)
Hi, please could anyone tell me how to get 2D hilbert transform spectrum plot (timenormalized frequency) ? I would be thankful
Francesco (view profile)
GARIMA SOOD (view profile)
hi i retried it and now in plotting hht i am getting the error that
"imf" previously appeared to be used as a function or command, conflicting with its use here as the name of a variable.
A possible cause of this error is that you forgot to initialize the variable, or you have initialized it implicitly using
load or eval.
can someone help me in resolving this?
GARIMA SOOD (view profile)
hi I am getting the error as the dimensions of the martices to be concatenated are not consistent .i am also getting an error in the spline function i fail to understand what is wrong with my code i'm giving the input x as wavread('Hum.wav'); and am using the order as :
finding peaks
emd
plot hht
Abdelilah Erraji (view profile)
Agata Szczepanska (view profile)
Hi, please could anyone tell me how to get 2D hilbert transform spectrum plot (timenormalized frequency) ? I would be thankful
Abdul Rahman Hassan (view profile)
How to enter the inputs (Data) (x&y) into the matlab so I can run the code.
The data are around 600 points (x&y) in CSV file.
Where x represents the time and y the time.
chunlin fu (view profile)
jing zhang (view profile)
Mikolaj Kegler (view profile)
Francois Brault (view profile)
I love the audio file in the package!!!
latebird (view profile)
super (view profile)
Rahi Shet (view profile)
u can use default function findpeaks as [p,n]=findpeaks(x) where n represents the indices of the peaks
Rahi Shet (view profile)
The builtin function findpeaks gives the peak values and not the indices of the peak values.
Adam Pink (view profile)
Is it possible to generate the finalised data set of the hilbert transformed IMFs (ie. the data that is generated and plotted in the plothht function)?
Zhangjun Yu (view profile)
The file 'emd.m' will call the file 'findpeaks.m', if I can delete it and use the builtin function 'findpeaks'?
Zhangjun Yu (view profile)
I used 'audioread' instead and it works well.
Zhangjun Yu (view profile)
I typed the example code "[x,Fs] = wavread('Hum.wav')" in the command window (Matlab 2016a), but it told me that Undefined function or variable 'wavread'. Is it a removed function?
Lauri Palokangas (view profile)
This library has proven greatly useful, thanks for providing it.
I have one question about the residue. It seems that whatever sample data I feed to the HHT function, the residue (the last plotted IMF) always seems to be a negative parabola. I have understood that the residue characterizes the overall trend of the data, so it appears odd that the shape is always the same regardless of the data.
Thanks in advance for your responses.
Elsaeed Ali (view profile)
JGGJ
Ketan R (view profile)
I'm newbie...
working on hht...
trying to run those files....
[
emd(x)
Empiricial Mode Decomposition (HilbertHuang Transform)
findpeaks(x)
Find peaks.
plot_hht(x,Ts)]
couldn't run ...
help me...
renu (view profile)
dear frnds
thanks for guiding me.
please help me how can i use wavelet transform for analsing ecg signals.
Tudor Sascau (view profile)
cheng ghee lim
cheng ghee lim (view profile)
I am having issue on how to process the hilbert huang spectral analysis. Below is my code for the emd, can somebody help me???
sig=csvread('1 bar.csv');
sig1=sig(:,1);%Time
sig2=sig(:,2);%Data
a=decimate(sig1,4);
b=decimate(sig2,4);
t=a;
x=b;
figure(1);
plot(t,x);
axis tight;
title('Original Response');
xlabel('time');
ylabel('amplitude');
%%%%%%EMD%%%%%%%%
imfs=eemd(x, 0, 1);
imfs=imfs'; % for emd and eemd only
num=size(imfs,1); % Number of IMFs
figure(2)
subplot(6,2,1);
plot(t,imfs(1,:));
axis tight
subplot(6,2,2);
plot(t,imfs(2,:));
axis tight
subplot(6,2,3);
plot(t,imfs(3,:));
axis tight
subplot(6,2,4);
plot(t,imfs(4,:));
axis tight
subplot(6,2,5);
plot(t,imfs(5,:));
axis tight
subplot(6,2,6);
plot(t,imfs(6,:));
axis tight
subplot(6,2,7);
plot(t,imfs(7,:));
axis tight
subplot(6,2,8);
plot(t,imfs(8,:));
axis tight
subplot(6,2,9);
plot(t,imfs(9,:));
axis tight
subplot(6,2,10);
plot(t,imfs(10,:));
axis tight
subplot(6,2,11);
plot(t,imfs(11,:));
axis tight
subplot(6,2,12);
plot(t,imfs(12,:));
axis tight
%%%%%%%Instantaneous frequency%%%%%%%
imfs=imfs';
omega_m1=ifndq(imfs(:,2),1);%IMF using for instantaneous frequency
omega_m2=ifndq(imfs(:,3),1);%IMF using for instantaneous frequency
omega_m3=ifndq(imfs(:,4),1);%IMF using for instantaneous frequency
omega_m4=ifndq(imfs(:,5),1);%IMF using for instantaneous frequency
omega_m5=ifndq(imfs(:,6),1);%IMF using for instantaneous frequency
omega_m6=ifndq(imfs(:,7),1);%IMF using for instantaneous frequency
omega_m7=ifndq(imfs(:,8),1);%IMF using for instantaneous frequency
omega_m8=ifndq(imfs(:,9),1);%IMF using for instantaneous frequency
omega_m9=ifndq(imfs(:,10),1);%IMF using for instantaneous frequency
omega_m10=ifndq(imfs(:,11),1);%IMF using for instantaneous frequency
omega_m11=ifndq(imfs(:,12),1);%IMF using for instantaneous frequency
omega_m12=ifndq(imfs(:,13),1);%IMF using for instantaneous frequency
figure (3)
subplot(4,2,1);
plot(t,imfs(:,2));
grid;
title('Intrinsic Mode Function 1');
subplot(4,2,2);
plot(t,imfs(:,3));
grid;
title('Intrinsic Mode Function 2');
subplot(4,2,3);
plot(t,imfs(:,4));
grid;
title('Intrinsic Mode Function 3');
subplot(4,2,4);
plot(t,imfs(:,5));
grid;
title('Intrinsic Mode Function 4');
subplot(4,2,5);
plot(t,imfs(:,6));
grid;
title('Intrinsic Mode Function 5');
subplot(4,2,6);
plot(t,imfs(:,7));
grid;
title('Intrinsic Mode Function 6');
figure (4)
subplot(4,2,1);
plot(t,imfs(:,8));
grid;
title('Intrinsic Mode Function 7');
subplot(4,2,2);
plot(t,imfs(:,9));
grid;
title('Intrinsic Mode Function 8');
subplot(4,2,3);
plot(t,imfs(:,10));
grid;
title('Intrinsic Mode Function 9');
subplot(4,2,4);
plot(t,imfs(:,11));
grid;
title('Intrinsic Mode Function 10');
subplot(4,2,5);
plot(t,imfs(:,12));
grid;
title('Intrinsic Mode Function 11');
subplot(4,2,6);
plot(t,imfs(:,13));
grid;
title('Intrinsic Mode Function 12');
figure(5)
plot(t, omega_m1/2/pi,'r',t, omega_m2/2/pi,'r',t, omega_m3/2/pi,'r',t, omega_m4/2/pi,'r',t, omega_m5/2/pi,'r',t, omega_m6/2/pi,'r',t, omega_m7/2/pi,'r',t, omega_m8/2/pi,'r',t, omega_m9/2/pi,'r',t, omega_m10/2/pi,'r',t, omega_m11/2/pi,'r',t, omega_m12/2/pi,'r');
grid;
title('All Instantaneous Frequency');
figure(6)
subplot(6,1,1);
plot(t, omega_m1/2/pi,'r');
grid;
title('Instantaneous Frequency IMF1');
subplot(6,1,2);
plot(t, omega_m2/2/pi,'r');
grid;
title('Instantaneous Frequency IMF2');
subplot(6,1,3);
plot(t, omega_m3/2/pi,'r');
grid;
title('Instantaneous Frequency IMF3');
subplot(6,1,4);
plot(t, omega_m4/2/pi,'r');
grid;
title('Instantaneous Frequency IMF4');
subplot(6,1,5);
plot(t, omega_m5/2/pi,'r');
grid;
title('Instantaneous Frequency IMF5');
subplot(6,1,6);
plot(t, omega_m6/2/pi,'r');
grid;
title('Instantaneous Frequency IMF6');
figure(7)
subplot(6,1,1);
plot(t, omega_m7/2/pi,'r');
grid;
title('Instantaneous Frequency IMF7');
subplot(6,1,2);
plot(t, omega_m8/2/pi,'r');
grid;
title('Instantaneous Frequency IMF8');
subplot(6,1,3);
plot(t, omega_m9/2/pi,'r');
grid;
title('Instantaneous Frequency IMF9');
subplot(6,1,4);
plot(t, omega_m10/2/pi,'r');
grid;
title('Instantaneous Frequency IMF10');
subplot(6,1,5);
plot(t, omega_m11/2/pi,'r');
grid;
title('Instantaneous Frequency IMF11');
subplot(6,1,6);
plot(t, omega_m12/2/pi,'r');
grid;
title('Instantaneous Frequency IMF12');
Dulip (view profile)
function hilbert is included in matlab signal processing toolbox. hope this helps.
Mhammad ali (view profile)
function hilbert is not defined.. can anyone explain??
nermine monla (view profile)
hi, please can anyone tell me how to get 3D plot (amplitudetimefrequency)after hilbert spectrum analysis of IMF?
omid (view profile)
Hum.wav
J.Hax (view profile)
I'm having issues with this using Matlab R2012a student version. The issue in the EMD file is the getspline command. Anyone else having the same issues with this?
Mahmoud Hassan (view profile)
renu (view profile)
I want to know how we can analyse a composite signal using HHT. and how HMS is used to analyse the signal. Please send me the code for finding Hilbert Marginal spectrum. I would be thankful to you
Baha (view profile)
Are intrinsic modes correspond to "natural modes" of a structure? For example, of a nonlinear time history of a building, if I use acceleration or displacement response and use this code HilbertHuang Transform, am I gonna get responses each corresponding to natural modes of the structure?
Thanks,
David Burton (view profile)
...and noise. Lots of noise.
I wrote, "there is probably an Nyear periodic component, some shorterperiod components, and a linear trend."
I should have written, "there is probably an Nyear periodic component, some shorterperiod components, a linear trend, and noise."
David Burton (view profile)
I have a question about the applicability of HHT/EMD.
Suppose I have a time series dataset in which there is probably an Nyear periodic component, some shorterperiod components, and a linear trend. Unfortunately, my dataset is considerably shorter than N years. My goal is not to find the periodic components (though that would be mildly interesting), but to more accurately find the linear trend, by removing the influence of the periodic components.
Is HHT/EMD useful for this? Can it help me make a more accurate determination of the linear trend than simple linear regression analysis would do?
Ching Kevin (view profile)
mira (view profile)
hi,
i have one question. what is the frequency label in timefrequency plot refer to? is it in Hz or rad/sec?
thank you.
lu (view profile)
Le (view profile)
kai (view profile)
Thanks!
dimple khanna (view profile)
lin (view profile)
Traian Preda (view profile)
Hi,
I have a question, regarding how I can identify in the last plot which frequencies are plotted in time?
Thank you
n d (view profile)
is this code helpful in video stabilization using HHT
francesco (view profile)
thank for the tool,
Alan Tan (view profile)
Hi Leon, thanks for pointing out. I tried findpeaks in the signal processing toolbox, and but it returned "No peaks found." Possibly, in the later versions (version I am using now is kinda old), the function has been corrected to locate the peak. Actually, I am not so concerned about the error because it only arises when there are instances of adjacent samples having exactly the same value on the left side of the peak, e.g., findpeaks([1 2 2.001 3 3 2 1]) returns the correct answer. In most cases, at least in the context of HHT/EMD or general signal processing, such occurrences are rare. Sorry I could not be more helpful.
Leon (view profile)
The "findpeaks" function is wrong. For signals such as [1,2,2,3,3,2,1], it return both 2nd and the 4th time points as local peaks, while only the 4th should be considered a peak. The MATLAB findpeaks function (in Signal Processing Toolbox), meanwhile, does not show this fallacy.
jianning (view profile)
I am a student in the university.Millions of thanks will go to you if you sent me this file.
junqua (view profile)
RAM KINKER MISHRA (view profile)
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 (view profile)
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 noncommercial, can HHT be used freely?
Please help me out.
Thanks
dedenlaut (view profile)
How to plot this result in energyfrequencytime domain?
Sergi Andruschenko (view profile)
RAM KINKER MISHRA (view profile)
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:end1)).^2)/norm(a(j,1:end1))^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 (view profile)
please provide me the matlab code of hilbert transform
RAM KINKER MISHRA (view profile)
My aim is to find the mean frequency derived via HilbertHuang 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 (view profile)
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:end1)).^2)/norm(a(j,1:end1))^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 (view profile)
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 (view profile)
hi please anyone tell how to apply this one on images esp of medical images.
plz..reply fast.
aami (view profile)
cud any one plz give me the link to the code for decomposing an image using EMD...
Thanx in advance...
Meenu
Lukasz (view profile)
Shouldn't you unwrap the angle th before sending it to diff?
Omar (view profile)
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 (view profile)
Bubo Bubic (view profile)
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 (view profile)
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 (slopebased 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 timefrequencyenergy distribution, you should plot abs(hilbert(imf{k})).^2 as the value on the Z axis (plot_hht.m line 15).
Alan Tan (view profile)
@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 (view profile)
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 (view profile)
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 timefreq 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 (view profile)
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 (view profile)
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 (view profile)
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 (view profile)
The following reply is given to Agnes over email.
Hi Agnes, The 2 initial plots are the timefrequency plots of the HHT corresponding to the 2 IMFs with of the highest energy. Subsequent plots are the IMFs, plotted against time.
Agnes (view profile)
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 (view profile)
Thanks for your sharing.
delinda he (view profile)
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 = 1b/max(b);
% Set timefrequency plots. I do not how it plot the timefrequency figures too
N = length(x);
c = linspace(0,(N2)*Ts,N1);
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 (view profile)
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 (view profile)
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 (view profile)
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 (view profile)
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 (view profile)
why can't i use it correctly?
The user would also need the Signal Processing Toolbox to perform the Hilbert Transform itself, though this is elementary to code.
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.
Forgot to add that you will need to have the Spline Toolbox to run this one.