Last time, I asked how I can use the command DiscreteWaveletPacketTransform[data, filter, 0] for spike detection and I got answers mostly based on an automatically defined threshold. In my case it is crucial to first define a threshold to distinguish between the real and spurious data, so built-in functions are not my favorites.

I developed my own version of code but as I use a filter like Haar, there would be no difference between a spike on point 19 or 20, as this filter works inherently with pairs of numbers. Once again, in Mathematica version 7 using the above command I had no problem! I read the links given already but no use. That would be great if someone can enhance my code such that it can differentiate between two consecutive spikes.

data = Table[Sin[x] + Random[], {x, 1, 10, 0.1}]
data[[20]] = 100; data[[40]] = 100;

enter image description here

dwt = DiscreteWaveletTransform[data, HaarWavelet[], 1];
Sigma = Mean[ Abs[dwt[[1, 2]] - Mean[dwt[[1, 2]]]]]/0.6745;
Lambda = N[Sqrt[2 Log[Length[dwt[[1, 2]]]]]];
thresh = Lambda*Sigma;
shr[c_, wind_] := If[Abs[#] >= thresh, 1, 0] & /@ c;
dwtS = WaveletMapIndexed[shr, dwt]
ListLinePlot[InverseWaveletTransform[dwtS, HaarWavelet[]]]

enter image description here


2 Answers 2


I don't know if you are wedded to wavelets but... Have you considered a "compound median filter" (q.v.) ?

For a list of data x and filter width 2r+1,

MedianFilterRoot[x_, r_] := FixedPoint[MedianFilter[#, r] &, x]  
CompoundMedianFilter[x_, r_] := 
 Fold[MedianFilterRoot[#1, #2] &, x, Range[r]]   

Plotting CompoundMedianFilter[x,r-1]-CompoundMedianFilter[x,r] shows all spikes of width r. Your example with two consecutive spikes could be something like:

data = Table[Sin[x] + Random[], {x, 1, 10, 0.1}]  
data[[20]] = 60; data[[21]] = 80; data[[40]] = 100;  

Running CompoundMedianFilter[data,0]-CompoundMedianFilter[data,1] returns two unit-width spikes of amplitude 20 and 100, at indices 21 and 40, respectively.
Similarly, CompoundMedianFilter[data,1]-CompoundMedianFilter[data,2] returns a width-2 spike of amplitude 59 at indices 20 and 21.
Hence, adjacent spikes are resolved with essentially undistorted amplitudes (for this example with relatively large spikes).

To reconstruct the spike-free signal, subtract all spikes of significant amplitude (and any width) from the original data.

  • $\begingroup$ Welcome to Mathematica.SE! Please consider registering your account so that any upvotes you get on this answer are added to those you might get on future questions and answers. That way, over time you will be able to do more on the site (post graphics, edit things, etc). $\endgroup$
    – Verbeia
    Commented Sep 1, 2012 at 8:03

Peak Detection Process may involve 3 main steps

  • smoothing
  • baseline correction
  • peak picking

For smoothing you can use: Moving average filter, Savitzky-Golay filter, Gaussian filter, Kaiser window, Continuous Wavelet Transform, Discrete Wavelet Transform, Undecimated Discrete Wavelet Transform

For baseline correction you can use: Monotone minimum, Linear interpolation, Loess, Continuous Wavelet Transform, Moving average of minima

Peak Finding Criterion you can use: Detection/Intensity threshold, Slopes of peaks, Local maximum, Shape ratio, Ridge lines, Model-based criterion, Peak width.

In your problem you need to do baseline correction in order to find distinctive peaks.

  • $\begingroup$ Smoothing is what you do to get rid of the spikes, not to detect them. $\endgroup$
    – stevenvh
    Commented Aug 26, 2012 at 15:13
  • $\begingroup$ Not really. The aim of smoothing is to give a general idea of relatively slow changes of value with little attention paid to the close matching of data values. check en.wikipedia.org/wiki/Smoothing $\endgroup$
    – s.s.o
    Commented Aug 26, 2012 at 18:23
  • $\begingroup$ Simplest smoothing filter is moving average. The following code decimates the peaks: data = Table[Sin[x] + Random[], {x, 1, 10, 0.1}]; data[[20]] = 100; data[[40]] = 100; llp1 = ListLinePlot[data, PlotRange -> All]; dataMA = MovingAverage[data, 10]; llp2 = ListLinePlot[dataMA, PlotRange -> All]; Show[llp1, llp2] $\endgroup$
    – stevenvh
    Commented Aug 26, 2012 at 18:50
  • $\begingroup$ As part of peak detection process - moving average can do smoothing - that's also what I've written. But, it does not decimates just averaged the point to near 10 data points. The peaks are still there and spread to those points. $\endgroup$
    – s.s.o
    Commented Aug 26, 2012 at 19:09
  • $\begingroup$ By decimating I mean reduce its amplitude dramatically. The moving average is just the simplest filter. How about this continuous time filter? It's just a first order smoothing filter, and it leaves just 1 % of the peak's amplitude. To detect a peak you need a high-pass filter, not the low-pass filter which a smoothing filter is. $\endgroup$
    – stevenvh
    Commented Aug 27, 2012 at 5:58

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.