12
$\begingroup$

I am trying to learn Hilbert-Huang transform (HHT) and its applications (especially in ECG signal processing). There are several HHT packages (scripts) HHT for R and Matlab. Is there any known solution for Mathematica or I should write my own implementation of HHT?

$\endgroup$
2
  • $\begingroup$ If there's an R implementation you could perhaps use RLink to make use of it. $\endgroup$ Jul 17, 2013 at 19:49
  • $\begingroup$ Or MATLink to call MATLAB from Mathematica. However, glancing through the MATLAB file it seems easy enough to translate it to Mathematica. $\endgroup$ Jul 17, 2013 at 20:43

1 Answer 1

19
$\begingroup$

Here a more or less straightforward translation of Alan Tan's MATLAB code to Mathematica code:

emd[x_List] :=
 Module[{xt = x, imf={}, x1, x2, s1, s2, sd},
  While[ ! isMonotonic[xt],
    x1 = xt;
    sd = Infinity;
   While[ (sd > 0.1) || ! isIMF[x1],
    s1 = getSpline[x1];
    s2 = -getSpline[-x1];
    x2 = x1 - (s1 + s2)/2;
    sd = Total[(x1 - x2)^2]/Total[x1^2];
    x1 = x2;
   ];          
   AppendTo[imf, x1];
   xt = xt - x1;
  ];
  Append[imf, xt]
]

isMonotonic[x_List] := Length[findPeaks[x]] Length[findPeaks[-x]] == 0;

isIMF[x_List] :=
 Module[{u1, u2},
  u1 = Count[Most[x] Rest[x], _?Negative];
  u2 = Length[findPeaks[x]] + Length[findPeaks[-x]];
  Abs[u1 - u2] <= 1
  ]

getSpline[x_List] :=
 Module[{n, p},
  n = Length[x];
  p = findPeaks[x];
  Interpolation[Transpose[{Flatten[{0, p, n + 1}], Flatten[{0, x[[p]], 0}]}]][Range[n]] 
  ]

findPeaks[x_List] :=
 Module[{n, u},
  n = Flatten@Position[Differences[Boole[# > 0] & /@ Differences[x]], _?(# < 0 &)];
  u = Flatten@Position[x[[n + 1]] - x[[n]], _?(# > 0 &)];
  n[[u]] + 1
 ]

Off[Interpolation::inhr];

Test:

SeedRandom[42];
testData = MovingAverage[RandomReal[{-10, 10}, {200}], 10];

ListLinePlot[emd[testData], PlotRange -> All]

Mathematica graphics

ListLinePlot[{testData, emd[testData][[4]], getSpline[testData], -getSpline[-testData]}]

Mathematica graphics

$\endgroup$
5
  • $\begingroup$ Wow. I really thought about doing the same, but the C variant. Impressing how much effort you've spent. (+1) $\endgroup$
    – Stefan
    Jul 18, 2013 at 7:50
  • $\begingroup$ @Sjoerd C. de Vries Thanks for doing this. (Sorry for two comments button problems.) I have been investigating. Unfortunately I have very long data sets and was wondering if you know of any method to speed up the calculation. I have tried changing the interpolation order (no effect). I have done some searching but not found any straightforward ideas. Any thoughts? Some ideas: i) split the data into short (overlapping?) lengths and apply to each length, ii) bandpass filter the data, iii) try to find the IMFs associated with the frequencies I am interested in. $\endgroup$
    – Hugh
    Dec 31, 2014 at 18:49
  • $\begingroup$ @hugh from the looks of the code I suppose some speed up can be achieved by getting rid of the Append and AppendTo which are known to be slow especially for large datasets. There are good alternatives for that. I feel Count could be done faster too, probably using he UnitStep trick. Unfortunately I don't have time to look into that right now given that it's New Years eve. $\endgroup$ Dec 31, 2014 at 20:35
  • $\begingroup$ @hugh I've analyzed this a bit further and it seems the main problem is the complexity of Alan Tan's algorithm. Counting the number of inner loop executions as a function of data list length it looks like it's O(N^4) with worst cases having even higher complexity O(x^4.5) or worse. Plot, Log plot. So, it seems there's not much one can do other than invent a better algorithm. $\endgroup$ Jan 1, 2015 at 14:33
  • 1
    $\begingroup$ @SjoerdC.deVries Thanks for looking at this. Not good news. I noted some speed up by improving the interpolation order but this is not a good way forward. I am thinking of smoothing my data first and hoping that will then reduce the complexity. Many thanks. Happy New Year. $\endgroup$
    – Hugh
    Jan 1, 2015 at 15:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.