Ytukyg

Suppose we have (an interval of) a time series of measurements:

We assume it can be explained as a "simple" underlying signal overlaid by noise. I'm interested in finding a good algorithm to estimate the value of the simple signal at a given point in time -- primarily for the purpose of displaying it to human users who know, more or less, how to interpret the signal but would be distracted by the noise.

For a human observer of the plot it looks very clearly like the underlying signal has a jump discontinuity at about $t=18$. But that's a problem for automatic noise removal, because the techniques I know are all predicated of a "nice" underlying signal meaning a "smooth" one. A typical anti-noise filter would be something like colvolving with a Gaussian kernel:

which completely fails to convey that the left slope of the dip is any different from the right one. My current solution is to use a "naive" rolling average (i.e. convolution with a square kernel):

whose benefit (aside from simplicity) is that at least the sharp bends in the signal estimate alert the viewer that something fishy is going on. But it takes quite some training for the viewer to know what fishy thing this pattern indicates. And it is still a tricky business to pinpoint when the abrupt change happened, which is sometimes important in my application.

The widths of the convolution kernels in the two examples above were chosen to give about the same smoothing of the pure noise (since I've cheated and actually constructed the sample data I'm showing as the sum of a crisp deliberate signal and some explicit noise). If we make them narrower, we can get the estimate to show that there's an abrupt change going on, but then they don't remove all of the noise:

I can't be the first person ever to face this problem. Does it ring a bell for anyone? I'd appreciate ideas, pointers to literature, search terms, a conventional name for the problem, whatever.

Miscellaneous remarks:

1. No, I cannot rigorously define what it is I want to optimize for. That would be too easy.




2. It will be nice if the smoothed signal can show clearly that there's no significant systematic change in the signal before the jump at $t=18$.




3. In the example I show here, the jump was much larger than the amplitude of the noise. That's not always the case in practice.




4. I don't need real-time behavior, so it is fine that the estimate of the signal at some time depends on later samples. In fact I'd prefer to find a solution that commutes with time inversions.




5. Basing a solution on outside knowledge about how the particular signal I'm looking at ought to behave is not an option. There are too many different measurements I want to apply it to, and often it's something we don't have a good prior expectation for.

I am working at a group where we study sudden drops in current (due to ion current blockades) in some particular measurement setup.

We detect these drops by calculating the moving average (and in fact, as of recently also a moving standard deviation, since the noise level varies over time in "bad" measurements) and selecting current data points whose drop is more than $5sigma$ (ie significant drop).

Any contiguous set of these points is called an "event", and to make sure we have the right kind of events (as opposed to noise), we further select "events" by integrating the difference of the event and the moving average, looking at its duration (ie. number of data points), and its amplitude, but I suppose you won't need all of those.

What you could do is detect events like we do, and then do interpolation or smoothing only between events (piecewise).

I must admit that we probably miss (don't detect) some "actual" events because of noise issues. Looking at your example, however, this should not be a problem for you and our code would detect it without hesitation.

One more word on our research because it's so much fun:

In fact the drops in current are so small, and the noise is so big, in our case, that we first do an 8 pole Bessel filter (typically at 10kHz) before we can even see these drops. But once the events have been detected, we can go back to the unfiltered signal and continue analysis there.

For an introduction to our research, see this paper.

If I understand you correctly, your current approach (Gaussian and moving average filters) is adequate for detecting overall changes (what akkkk calls "events") but not for displaying the shape of the underlying data. Obviously, you can't truly know the shape of your overall data because it's masked by noise, but the eye seems to see it, so, the question is, shouldn't we be able to automatically show more or less the same line that the eye sees.

There are some techniques for this, but you'll have to experiment to see what works. I'm sure there are other techniques, but here's what I can think of:

• Broad-band noise reduction is used in audio to remove noise. It works well if your noise is broad spectrum and your signal is narrow spectrum, which seems to be the case here. It is tricky to implement. I don't have any good references for you, but maybe this book. I believe it's also covered in DAFX, but there have got to be better sources on it.




• Median filtering is an excellent technique that works well in certain kinds of image processing and would probably work well here. The mathematics behind it don't work out well (ie it's hard to prove it works), but in practical situations it often works well, and that's good enough for many people. I suspect it will work well here, but you'll have to try it.




• The gaussian filter and moving average filter you've used can be though of special cases of low-pass filters. If you have some idea of the frequency of your signal vs the frequency content of your noise (or if you are willing to experiment with your real data until your output corresponds to your eyes) you can design a low-pass filter. This is what akkkk was suggesting when he talked about the 8 pole Bessel filter (I'm not sure this particular filter would preserve the shape of the signal, but it probably would). Even a second order recursive filter would be a lot better than these filters, I imagine. I have a tutorial on simple audio filters to get you started.




• There exist adaptive filters which, depending on your situation, you might be able to use. See Wiener and Kalman filtering. I've never used a Kalman filter, but skimming the Wikipedia page it looks like it might be purpose-built for the kind of thing you are doing.

If these techniques seem daunting to you, try median filtering first. It's simple to implement and will probably help significantly.

Have two algorithms: one that checks for discontinuities, the other that takes each "smooth" section and uses a typical noise removal algorithm. For the first algorithm, I would recur through the entire time interval, first looking at very large chunks but getting smaller and smaller with each recursion, comparing the slope of each chunk's linear fit.