I am a new user of Mathematica (ver. 9.0) and I kindly ask your help for some data analysis. I have a list of measurement (Voltage) acquired at a constant frequency. This signal is interrupted periodically by a characteristic drop, which is intrinsic in the measurement method (see example below). The time zero is set at the middle of the first "drop". enter image description here

Every drop is constitued by several datapoints, i.e., is not a single measurement (second image). I would like to exclude the drops from the data to work on the real signal which is the one comprised between each fall in the intensity. Can you suggest an easy way to filter or delete this artificial signal variation?

Thanks a lot! enter image description here

P.S. I also provide a data sample, which you can download from this link: https://www.dropbox.com/s/0qspq669ksk7lai/Iamorph.txt?dl=0

  • 1
    $\begingroup$ Try MedianFilter[myData, 100]. $\endgroup$ – David G. Stork Apr 7 '19 at 23:10
  • 3
    $\begingroup$ Can you provide sample data to play with? (Also just out of curiosity, this looks like drops from opening the cover on a fluorimeter to add reagents?) $\endgroup$ – MelaGo Apr 7 '19 at 23:12
  • $\begingroup$ Do you want to remove the drops AND smooth the data? Or just remove the drops? It looks like once the drops to zero are removed there still exist jumps from one level to the next. Are those jumps to be smoothed or left alone? $\endgroup$ – JimB Apr 7 '19 at 23:38
  • 2
    $\begingroup$ OK. Then you definitely don't want to use MedianFilter. And to second @MelaGo 's suggestion: providing some sample data would get you quicker and better answers. $\endgroup$ – JimB Apr 8 '19 at 4:53
  • 1
    $\begingroup$ @MelaGo That was a nice guess! Actually the data come from trying to do spectroscopy during 3D printing $\endgroup$ – Dario Cavallo Apr 8 '19 at 5:03

This approach is based on the standard deviation of a moving window. It doesn't require regular periodicity of the drops, but the standard deviation cutoff and moving-window size may need to be adjusted.

sdcutoff = .1;
movingwindow = 8;
mw := If[EvenQ[movingwindow], movingwindow, movingwindow + 1] ;
sd = Table[{data[[n, 1]], 
     data[[All, 2]][[n - mw/2 + 1 ;; n + mw/2]]]}, {n, mw/2, 
    Length[data] - mw/2}];

Plotting the standard deviation shows the spikes that coincide with the drops (and helps select an appropriate sdcutoff).

ListPlot[sd, PlotRange -> All, AxesLabel -> {"Time [sec]", "Standard Deviation"}]

enter image description here

Only keep data[[t,y]] for which the standard deviation of the moving window centered at t is < sdcutoff:

dataculled = {};
Do[If[sd[[n - mw/2 + 1, 2]] < sdcutoff, 
   AppendTo[dataculled, data[[n]]]],
  {n, mw/2, Length[data] - mw/2}];

ListPlot[{data, dataculled}, PlotStyle -> {Blue, Red}]

enter image description here

Here is another example of data with similar drops, but not at regular intervals (from a fluorescence anisotropy binding experiment):

ListPlot[data2, Joined -> True, AxesOrigin -> {0, -.1}]

enter image description here

sdcutoff = .03;
movingwindow = 30;
sd2 = Table[{data2[[n, 1]], 
     data2[[All, 2]][[n - mw/2 + 1 ;; n + mw/2]]]}, {n, mw/2, 
    Length[data2] - mw/2}];
ListPlot[sd2, PlotRange -> All]

enter image description here

dataculled2 = {};
Do[If[sd2[[n - mw/2 + 1, 2]] < sdcutoff, 
   AppendTo[dataculled2, data2[[n]]]],
  {n, mw/2, Length[data2] - mw/2}];
ListPlot[{data2, dataculled2}, PlotStyle -> {Blue, Red}]

enter image description here

| improve this answer | |
  • $\begingroup$ +1 This must be much more robust by not having to assume a specific periodicity. $\endgroup$ – JimB Apr 9 '19 at 2:56
  • $\begingroup$ @JimB The tradeoff is that you have to fool around with the sdcutoff and window size, and as you can see, there are still a few erroneous points included. There's probably a way to automate it more, maybe by iterating and looking at the variance of the sd's,.. but I would have to take more time to think about it. $\endgroup$ – MelaGo Apr 9 '19 at 3:14
  • $\begingroup$ I guess the approach is dictated by what one really knows about the process in addition to the data. If it's the "drops" and the drops have some consistent shape (as opposed to a consistent periodicity), then that knowledge could be exploited. In short, it requires more (subject matter) information than what the data contains. And many times getting that kind of information is hard to extract from whoever is asking the question. $\endgroup$ – JimB Apr 9 '19 at 3:57
  • $\begingroup$ @JimB, true. I was motivated to think about a solution without periodicity, because I sometimes have similar-looking data, such as in the second example. In that case, the "drops" (which can also be spikes) are due to a person opening the lid on the fluorimeter to pipette something into the cuvette. So obviously the interval and duration of the artifacts can vary quite a lot. $\endgroup$ – MelaGo Apr 9 '19 at 5:12

Let's define the scan-window from the periodicity of your data in the array data:


Let's scan your data with such window width:

ress = {};
 dd = data[[window*i + 1 ;; window*(i + 1)]];
 mm = MinMax@dd[[60 ;; 120, 2]];
 delta = mm[[1]]/mm[[2]];
 ress = Join[ress, Cases[dd, s_List /; s[[2]] > delta*mm[[2]]]],
 {i, 0, Floor@(Length@data/window) - 1}]

The result looks like this:

ListPlot[{data, ress}]

enter image description here

There are some mistakes of selector possible but it can be adjusted by the right choice of window size. It should be a bit bigger than the period of your data. An additional parameter is delta. Playing with it, you can change the selection results too.

| improve this answer | |
  • $\begingroup$ Thank you very much for your help! I will soon apply your solution! All the bests $\endgroup$ – Dario Cavallo Apr 8 '19 at 14:13

This method, too, relies on assuming a drop approximately every 218 points. (Maybe something not dependent on assuming a constant cycle is needed. If so, please indicate that in your question.)

(* Find the potential minimums approximately every 218 points and \
within 22 points of either side of the initial potential minimum *)
d = 22;
low = Flatten[
   Table[Select[data[[Range[Max[1, i - d], Min[Length[data], i + d]]]],
     #[[2]] == 
       Min[data[[Range[Max[1, i - d], Min[Length[data], i + d]], 
          2]]] &],
    {i, 1, Length[data], 218}], 1];

(* Get a list of indices d points to the left and d points to the \
right of each potential minimum *)
pos = Flatten[
  Position[data, #] & /@ low]; (* Position of potential minimums *)
indices = 
    Table[Range[Max[1, pos[[i]] - d], Min[Length[data], pos[[i]] + d]],
     {i, Length[pos]}]]]];

(* Keep just the "good" data *)
goodData = Delete[data, Table[{indices[[i]]}, {i, Length[indices]}]];

(* Show results *)
Show[ListPlot[goodData], ListPlot[data[[indices]], PlotStyle -> Red]]


| improve this answer | |
  • $\begingroup$ Great!! Thanks a lot for your help, I will try to implement it ASAP! Bests $\endgroup$ – Dario Cavallo Apr 8 '19 at 14:12

Another approach:

data = Import["Iamorph.txt", "Table"];

(* select elements for a max distance to nearest neighbor *)
neigbors[limit_] := 
 Select[data, EuclideanDistance[#, Nearest[data, #, 2][[2]]] < limit &]

towns = neigbors[.02];

ListPlot[{data, towns}, 
 PlotStyle -> {Black, Directive[PointSize[.005], Red]}]

enter image description here

(* eliminate the sprinklink of small values *)
bigTowns = Select[towns, #[[2]] > .1 &]; bigTowns // ListPlot

enter image description here

(* moving median eliminates most outliers *)
medians = MovingMedian[bigTowns, 5]; medians // ListPlot

enter image description here

(* find the towns *)
clusters = 
  FindClusters[medians, 11, DistanceFunction -> EuclideanDistance, 
   Method -> "Agglomerate"];

enter image description here

(* and plot the mean values for each *)
Mean /@ clusters // ListPlot

enter image description here

| improve this answer | |

Below I am reusing the approach / code from this answer.


data = Import["~/Downloads/Iamorph.txt", "TSV"];

cleanData =
  Fold[Function[{data, pars},
       QRMonEcho[Style[pars, Purple, Bold]]⟹
       QRMonQuantileRegression[pars["knots"], pars["fraction"]]⟹
       QRMonErrorPlots["RelativeErrors" -> False]⟹
    ], data, {<|"knots" -> 50, "fraction" -> 0.5, "threshold" -> 1|>, <|"knots" -> 50, "fraction" -> 0.5, "threshold" -> 0.3|>}];

enter image description here

enter image description here

ListPlot[cleanData, PlotTheme -> "Detailed"]

enter image description here

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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