16

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

Community
  • 1
K-1
  • 581
  • 1
  • 3
  • 12

2 Answers2

4

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.

withparadox2
  • 2,481
  • 1
  • 24
  • 28
KennyColnago
  • 15,209
  • 26
  • 62
  • 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). – Verbeia Sep 01 '12 at 08:03
2

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.

s.s.o
  • 4,559
  • 2
  • 27
  • 42
  • Smoothing is what you do to get rid of the spikes, not to detect them. – stevenvh Aug 26 '12 at 15:13
  • 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 http://en.wikipedia.org/wiki/Smoothing – s.s.o Aug 26 '12 at 18:23
  • 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] – stevenvh Aug 26 '12 at 18:50
  • 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. – s.s.o Aug 26 '12 at 19:09
  • 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. – stevenvh Aug 27 '12 at 05:58