Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
If there's an R implementation you could perhaps use RLink to make use of it. –  Sjoerd C. de Vries Jul 17 '13 at 19:49
    
Or MATLink to call MATLAB from Mathematica. However, glancing through the MATLAB file it seems easy enough to translate it to Mathematica. –  Sjoerd C. de Vries Jul 17 '13 at 20:43

1 Answer 1

up vote 11 down vote accepted

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

share|improve this answer
    
Wow. I really thought about doing the same, but the C variant. Impressing how much effort you've spent. (+1) –  Stefan Jul 18 '13 at 7:50
    
@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. –  Hugh Dec 31 '14 at 18:49
    
@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. –  Sjoerd C. de Vries Dec 31 '14 at 20:35
    
@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. –  Sjoerd C. de Vries Jan 1 at 14:33
    
@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. –  Hugh Jan 1 at 15:51

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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