License

Copyright (c) 2016, Alan Tan
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.

Comments and Ratings (83)

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.

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

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.

%%%%%%%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

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 Hilbert-Huang Transform, am I gonna get responses each corresponding to natural modes of the structure?

Suppose I have a time series dataset in which there is probably an N-year periodic component, some shorter-period 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?

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.

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.

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

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.

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

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

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.

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

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

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?

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?

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.

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!

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.

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

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.

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