Need a better step detection algorithm

Question

I have a time series with lots of steps/jumps (data file here). A plot is given below. I would like to subtract an appropriate value for each of these square wave features to bring them back down to the baseline of the signal. A median filter works really well for removing a small number of outliers in a row, but in this case I probably need a different approach since the square wave jumps can have different durations as seen. A common method I've seen for doing this is to compute first differences of adjacent samples, and look for large differences to detect jumps. I implemented this method but the problem is it often fails, since the one tunable parameter for the method is a threshold value $t$ which the first differences must cross in order to detect a jump: $$ | x_{i+1} - x_i | > t $$ As can be seen in the plot below, the jumps I have are often different sizes, so a constant threshold value isn't the best approach. In particular, in some cases there is an interesting signal where adjacent samples can change by large values without being a jump! I have highlighted such a region in red.

Below is a zoomed in view of the red box area. You can see there is a square wave jump followed by an interesting signal. The red arrows depict a place where adjacent samples from an interesting signal have a larger distance between them than some of the jumps in the signal. Therefore a constant threshold method with finite differencing will not work for me.

Does anyone know of a robust procedure to detect and subtract the square wave jumps to end up with a smoothly varying signal with no jumps? I'm sure this must be a solved problem but I haven't had much luck searching online.

Answer 1

This is an interesting problem. I've downloaded your data and written a small program to process it.

There are missing data points that don't seem to be correlated with the jumps.

It seems to me the best approach is to look for matching jump values and adjust the values in between. The first step was to measure the RMS of the jumps from point to point to gauge the typical jump size. Next, I identified the jumps that exceed four times the RMS.

These are the results:

The Jump RMS = 1.270752

1262310195, -29.989300 Down
1262310259, 30.380400 Up
1262310355, -7.228400 Down
1262310364, 5.666300 Up
1262310438, 6.014400 Up
1262311415, 30.857000 Up
1262311626, -14.630000 Down
1262311627, -15.643500 Down
1262312656, -16.754800 Down
1262312657, -10.732500 Down
1262312859, 29.599500 Up
1262313801, 5.085400 Up
1262314529, 31.463000 Up
1262314694, -30.519800 Down
1262316412, -29.813200 Down
1262316483, 29.996800 Up
1262316604, -5.110200 Down
1262316612, -5.084100 Down
1262316614, 8.913200 Up
1262316618, -10.073700 Down
1262316620, 10.929100 Up
1262316621, 6.930200 Up
1262316623, -5.136700 Down
1262316624, 8.932800 Up
1262316625, -6.066600 Down
1262316628, 17.934800 Up
1262316630, -7.129600 Down
1262316631, -11.064500 Down
1262316659, 6.889600 Up
1262316661, -5.111200 Down
1262316663, -8.125000 Down
1262316665, 11.254100 Up
1262316667, -9.666300 Down
1262316668, 6.854200 Up
1262316680, -7.209700 Down
1262316689, -7.270200 Down
1262317545, 29.889600 Up
1262317785, -29.615400 Down
1262318838, -17.595400 Down
1262318839, -12.571100 Down
1262319054, 29.935400 Up
1262319539, -17.311700 Down
1262319540, -11.272300 Down
1262319567, 30.747600 Up
1262320679, 30.134800 Up
1262320869, -30.649600 Down
1262322664, -29.883400 Down
1262322703, 30.569600 Up
1262323699, 29.895000 Up
1262323975, -29.298600 Down
1262325018, -17.261100 Down
1262325019, -11.246300 Down
1262325269, 15.262600 Up
1262325270, 13.273800 Up
1262326865, 30.136500 Up
1262327071, -29.004400 Down

It looks like now it is just a matter of matching up and down swings. There is a wrinkle in that it appears that some of the jumps in the same direction are adjacent, so it takes two data points to make the jump.

Another wrinkle is when the jumps are matched, the values are quite exactly the same, so coming up with a formula for the proper rectifying shift is another puzzle.

Anyway, send me an email to cedron at protonmail dot com and I will send you the source code so far, and maybe take further yet.

Ced

Answer 2

This is not exactly an answer, but an attempt to illustrate further the problems I'm having with the first differences method. In the figure below, I focus on that same region as in the original question. I have added a median filtered curve in green (7-point centered window), plus the curve of first differences of the median filtered curve in red ($x_{i+1}-x_i$). Finally in blue I plot a constant "threshold" value of $\pm 4$ to identify a step.

You can see that the first difference curve (red) identifies the step jump quite well, but around 30:00 the red curve also passes the threshold value when in fact there is no jump there. I cannot increase the threshold value any more because of the following situation a little later in the time series: Here you can see two step features in a row whose first differences don't quite pass the $\pm 4$ threshold. So these jumps are not correctly identified. So the question becomes how to "enhance" the first difference signal in the presence of a step, while reducing it in non-step cases like in the first figure above. One idea I have is to compute the median absolute deviation (MAD) of the previous N samples, and compare the first difference value with the MAD instead of using a constant threshold. But I haven't gotten this method to work very well either, because sometimes the MAD is very small, and even a modest first difference value would result in false positive.

An updated data file with the median filter and first difference values is here. I would greatly appreciate any ideas!

Answer 3

I managed to get a fairly reliable solution to this problem using the following steps:

1. Smooth and differentiate the signal with a 2nd order Gaussian filter: $$ y(t) = \frac{d^2}{dt^2}G_{\sigma}(t) * x(t) $$
2. Search for zero crossings in $y(t)$ and record their positions (a first derivative filter will show peaks at each edge/jump. A second order filter will have zero crossings at the positions of the edge/jump).
3. Loop through the detected zero crossings $i$ (i.e. edge positions) and push the magnitude $|x_{i+1}-x_i|$ into a double-ended queue. I keep a cumulative sum of all elements in the queue extending backwards in time by some maximum amount (dt_max). If the next jump magnitude $|x_{i+1}-x_i|$ is within some relative error tolerance of the cumulative sum, the algorithm determines that a complete sequence of jumps has occurred leading back to the original signal baseline.
4. When a complete jump sequence is found, offsets are calculated for each portion of the jump to bring all the segments down to the signal baseline. They are popped off the back of the double-ended queue as they are processed.

This procedure is not perfect, and still misses a few jumps, but manages to capture most of them, even in the presence of significant noise.