|
I have a stream of data coming from a respiration belt and I need a good, robust algorithm to help me identify the peaks and troughs as the data is coming in (i.e. live data, not later analysis). The data is roughly sinusoidal with a period of a maybe 2 - 4 seconds. The data is polled about once every 100 ms. The data is noisy, but it's not too bad. We sometimes see sudden spikes if the person wearing the belt fidgets around, but if they are sitting still, it's fairly clean.
Obviously you can't identify a peak or a trough until after it's happened, but I'm trying to find the best and fastest way to identify them so we can use them as a trigger for something else. I understand that there are trade-offs involved.
I tried calculating an instantaneous derivative by comparing the change in signal between the current point and the second to last point and then looking for when that value falls within some threshold. That kinda works, but is a bit touchy. I tried to improve it by checking that when a point falls within the threshold it also has the opposite sign to the last identified peak/trough to try and eliminate having multiple points right around the peak identified as the peak. That helps a little. I also tried introducing, basically a filter, such that when a peak (or trough) is identified, I stop looking for another one for some period of time (say 1 second).
Does anybody have any better ideas?
|
|
|
|
|
Wjousts wrote: I stop looking for another one for some period of time (say 1 second).
Rather than that, you could wait for the slope (aka, instantaneous derivative) to reach a high value before allowing another low value to be detected. You can call this wait mode.
Wjousts wrote: I tried calculating an instantaneous derivative by comparing the change in signal between the current point and the second to last point and then looking for when that value falls within some threshold.
I'm not sure what you mean by "second to last point", but I'd say rather than looking for as instantaneous of a derivative, you look for an average slope of a few points. For example, take 4 data points to calculate a slope. Take a few of those in a row to ensure the correct pattern is being followed. So, if the current slope is less than the previous slope and the previous slope is less than the slope previous to that, then you can be reasonably certain that slope is actually decreasing. Basically, you want to be sure the data is trending as you expect it to be. When it falls below a threshold, then your minimum has been reached. Also, make sure to check if the slope has passed the minumum (e.g., with a negative slope followed by a positive one). Once that happens, you can enter wait mode before detecting a low slope again.
Wjousts wrote: Does anybody have any better ideas?
One more idea. Log some real data and save that so you can perform multiple tests on it. Perhaps start by graphing it out. Look for deviations from the expected pattern so you can account for them (e.g., the spikes you mentioned). Keep running your algorithm through that sample data until you have it working.
|
|
|
|
|
AspDotNetDev wrote: I'm not sure what you mean by "second to last point"
What I'm trying to do is basically this[^].
So, given a new point, the nth point, I can calculate the gradient at the previous point, n-1, but dividing the difference between the value of n-2 (the second to last point) and n with the difference in time between n-2 and n. I.e. I get an approximation of dV/dt. According to Wikipedia, the three point approximation is somehow better than the two point version?
Actually looking at this again made me realize that I've been attributing the gradient to point n, when in fact it belongs to n-1.
|
|
|
|
|
as it is more or less periodic, consider a Fourier analysis (google FFT). Keep the most recent 2^n samples, watch for the fundamental frequency and especially its phase.
Warning: if the harmonics are sufficiently large, the peaks and troughs will not coincide with the fundamental's phase, but then, which is it you actually want?
Luc Pattyn [My Articles] Nil Volentibus Arduum
The quality and detail of your question reflects on the effectiveness of the help you are likely to get. Please use <PRE> tags for code snippets, they improve readability. CP Vanity has been updated to V2.3
|
|
|
|
|
I've seen FFT before, but I didn't think it was useful for what I'm trying to do. I always thought it was more useful for post hoc analysis of data. I need to decide, as close as possible to when the event happens, when I hit a peak (and/or trough). Is there a way to use FFT for that? If I took a window of a few points (maybe 5 - 10), can I do something with FFT to tell me if there is a peak or trough in there?
|
|
|
|
|
If your signal is rather sinusoidal, I would definitely go for Fourier analysis. Take a sliding window, say 64 samples (that is 6.4 seconds, anywhere between 1 and 3 periods I guess); I would not go for a period shorter than 1 period. Calculate FFT, note the base frequency and phase (and maybe check the harmonics are small indeed). After 100 msec, drop the oldest sample, add the new one, and repeat.
If your signal is not sinusoidal at all, and you need the real peak/trough, then the above would not be good enough; I'm not sure anything else would under those conditions though.
Do you have some graphs to show us? And how accurate are the time points? is the sample period close to 100 msec, or is there a lot of jitter (say 100 msec plus/minus 20 msec). More jitter would cause very reduced accuracy. If there is significant jitter and the actual time is known, there may be ways to take advantage of that too.
Luc Pattyn [My Articles] Nil Volentibus Arduum
The quality and detail of your question reflects on the effectiveness of the help you are likely to get. Please use <PRE> tags for code snippets, they improve readability. CP Vanity has been updated to V2.3
|
|
|
|
|
Here's a graph showing a slice of some fairly clean data (look at the purple line - the blue line is the three point derivative, the triangles are one attempt at peak picking which actually doesn't look too bad ):
http://i55.tinypic.com/2uhwiyv.png[^]
I'd post the source data, but I don't have a good place to put it (suggestions are welcome).
The time points seem to be fairly accurately spaced at about 109 ms, but it relies on the accuracy of comparing DateTime objects in .NET, so it's plus or minus whatever that resolution is.
|
|
|
|
|
Ah, it is good to see some actual data. What isn't clear is this:
- what is the horizontal scale, what goes from 6 to 34? is this seconds?
- did you fit a curve? or is it just consecutive input points that got connected using straight lines?
Either way, it seems quite well behaved. The frequency anomaly in [16, 18] and the level dropping in [6, 18] would not be dealt with well when using Fourier analysis.
Now what is the reaction time you need? how many msec are you willing to wait before the latest peak/trough gets signaled?
(see another part of this thread). If a 300 msec delay is no problem, I would definitely go for a convolution filter.
Luc Pattyn [My Articles] Nil Volentibus Arduum
The quality and detail of your question reflects on the effectiveness of the help you are likely to get. Please use <PRE> tags for code snippets, they improve readability. CP Vanity has been updated to V2.3
|
|
|
|
|
Sorry, I was too lazy to label the axes. I used to beat up on students over that when I was in grad school.
Yes the x-axis is time in seconds. I cut the first 6 seconds off because it was a bit crappy (for some reason the signal was saturated for a few seconds at the beginning for reasons unknown). The vertical scale is arbitrary. It's proportional to voltage after it's passed through an ADC to give me a value between 0 - 4096.
I think 300 ms would probably be acceptable (but much more would probably be pushing it). Basically I'd like to do something during exhalation so that it'll be ready by the next inhalation. I played with the convolution filter a little and it does help remove some of the smaller shoulders or spurious peaks, so I think I might use it in combination with some other ideas and see what I end up with. If it helps, I think false positives (i.e. wrongly IDing a peak when it isn't one) is probably worse than false negatives (i.e. missing a real peak) so I think I might fiddle with it to see if there is anywhere I can trade off those objectives.
Edit: forgot to add, no the line isn't fitted, it's just the points connected up by Excel.
|
|
|
|
|
OK, this is what I would try:
1.
first use a kernel size 5, symmetric, sum of coefficients being 0 (i.e. some negative numbers, and a large positive center); example would be: (1, -6, 10, -6, 1). On a straight line this yields all zeroes no matter the slope. If not satisfied, fiddle with the kernel values; if that doesn't work, go for a larger kernel, which means more calculations (no problem), and more delay (could become a problem). For wider kernels, if you would need them, one often approximates a sin(x)/x function (i.e. a decaying periodic signal).
2.
to avoid false positives, I would postprocess with a low-pass, i.e. once you accept a peak/through at time T, reject anything that comes within the time period [T, T+delta] where delta has to be chosen carefully, based on what you consider a normal signal.
Luc Pattyn [My Articles] Nil Volentibus Arduum
The quality and detail of your question reflects on the effectiveness of the help you are likely to get. Please use <PRE> tags for code snippets, they improve readability. CP Vanity has been updated to V2.3
|
|
|
|
|
I've had some experience in processing moderately-to-EXTREMELY noisy data, and one technique that works for peak-finding is what the signal processing fraternity call a matched filter. Essentially, it's a block of coefficients that you walk along your data doing a point-by-point multiply-and-add, then apply a threshold (fixed or adaptive). To find a peak about 3 samples wide, your filter might look like [+1 +3 +5 +3 +1]. If you don't want another one too near, just tack some negative values on each end... There are all sorts of ways of theoretically developing filters, but I've had success with a pragmatic approach - get a good sample of clean data and match that.
Cheers,
Peter
Software rusts. Simon Stephenson, ca 1994.
|
|
|
|
|
convolution filters like the one you suggested are fine for a lot of applications including image processing, however I'm not sure they qualify as real-time: for a kernel size of 5, it takes 2 more samples for the peak to get into the middle of the kernel, and then maybe one more for the filter output to decrease again, signalling a peak was had. So that would cause a 200 to 300 msec delay in all. It may or may not be acceptable; if it is, it clearly would be the simplest approach.
Luc Pattyn [My Articles] Nil Volentibus Arduum
The quality and detail of your question reflects on the effectiveness of the help you are likely to get. Please use <PRE> tags for code snippets, they improve readability. CP Vanity has been updated to V2.3 modified on Wednesday, May 25, 2011 8:59 PM
|
|
|
|
|
Thanks for the word 'convolution'. Brain wasn't in gear - too much gardening. 
Real-timeness is, as you note, relative. I can't see an FFT beating a convolution filter for timeliness (given reasonable sizes for each).
Cheers,
Peter
Software rusts. Simon Stephenson, ca 1994.
|
|
|
|
|
Thanks. Interesting idea. I tried it on a slice of real data and it does smooth out some of the noise quite nicely. The problem still remains with picking up the actual peaks though. I guess I neglected to mention in my original post that the reason I can't just do something cheap and simple like using a threshold (or two) is because the baseline has a nasty habit of drifting up and down. I previously had an entirely different approach to this problem which used a flow sensor to measure the air being blown out by somebody as they breath and I used two thresholds there (a rising and a falling threshold to introduce hysteresis and cheaply filter out a good amount of low amplitude noise). It worked okay, but there where technical problems with reliably positioning the flow sensor.
|
|
|
|
|
re: Your baseline drift problem.
That's what I was referring to when I mentioned adaptive threshold. There are two simple ideas that spring to mind.
(If this reads like 'thinking aloud', it's meant to. Showing the thought process.)
1. Keep a moving average, which should be a reasonable measure of the current baseline
moving_average = weight * this_sample + (1.0 - weight) * moving_average
The value of weight adjusts the 'time constant' of the averaging, in units of the sample interval. In your case, a value in the range 0.01 to 0.1 is probably a good starting point. Then your threshold is moving_average + (an educated guess of say 75% of a typical peak)
2. This last bit looks messy, so why not run a moving peak detector ('leaky' peak detector) as well as a moving average as above, then you could get a good estimate of the current behaviour of your input signal. Both the moving average and leaky peak are very cheap to implement (in terms of processing effort).
Leaky peak:
if (this_sample > running_peak)
running_peak = this_sample;
else
running_peak = alpha * moving_average + (1.0 - alpha) * running_peak;
alpha probably about half of weight from above.
Do all this after smoothing, of course. Otherwise your noise will stuff up the running peak.
Cheers,
Peter
[edit] Fixed up leaky peak stuff. [/edit]
Software rusts. Simon Stephenson, ca 1994.
modified on Monday, May 30, 2011 7:11 PM
|
|
|
|
|
Thanks Peter. That looks pretty interesting. I played around with my one set of test data and it looks like it might work. I have to play around with the constants a little to see if I can get optimal values.
|
|
|
|
|
You were on the right track with instantaneous derivative.
At the actual peak the derivative goes to zero, one side of it is positive, the other side is negative. Remember though that your sample may not be at the peak, so the peak may be between 2 samples rather than the middle sample of a group of 3.
So all you need do to is find where the derivatives changes sign, pos to neg is peak, neg to pos is trough. It doesn't tell you the exact point of the peak/trough, just that one happened between the sign changes.
Hope this is some use
Steven
|
|
|
|
|
Thanks for the input Steven. That was pretty much what I started out doing, but it provided to be less robust than I wanted with the noise in the data. Some filtering (based on several other suggestion in this thread) seemed to help and then this method works a lot better.
|
|
|
|
|
Wjousts,
your signal indeed looks nice.
If it was perfect, all you would have to do is indeed look for changes in the sign of the derivative, i.e. detect monotonous sequences (increasing then decreasing, in alternation).
Added noise causes local perturbations of this dream situation, resulting in small breaches of the monotonocity. The closer you get to the peak, the more likely they get (as the derivative gets smaller and smaller).
My way to deal with that is to consider "quasi-increasing sequences" (resp. decreasing), i.e. values that go increasing but allow a backtracking limited by a threshold value. For example, assuming a threshold of 10, the following sequence is quasi-increasing: 0, 22, 31, 27, 45, 63..., while this one is not: 0, 22, 31, 20, 45, 63... (said differently, it stops being quasi-increasing at the value 20).
This approach is more robust than mere derivative computation and you can also filter on the total increase (decrease) of the sequence. You can also limit the length of the allowed backtracking to a given number of samples.
When will it detect a peak ? When the signal value has decreased to the peak value minus the threshold.
I guess that in any case you cannot avoid having to wait some times after a peak, to be sure it is a true one.
|
|
|
|
|
Thanks YDaoust. I like that idea a lot. It looks like it should be able to deal with shoulders or split peaks fairly nicely so long as a sensible threshold is used. I'll experiment with that idea a little.
|
|
|
|
|
I think so too. Your signal is clean, you should get good results. The nicety with the approach is that you work on the signal, not its derivative.
The threshold value can be set to a small multiple of the noise amplitude.
|
|
|
|
|
Removing spikes could be easily done by an FFT (converting to frequency spectrum), a low pass (removing high frequency components) and a reverse FFT 8going back to the time signal again). However, this pretty much kills your response time.
There are low-pass filter implementaitons in the time domain (wikipedia 1[^], wikipedia 2[^])
(I don't know wnough to recommend a specific algorithm, I am just parroting what "the guys" do.)
|
|
|
|
|
Thanks Peter. I've used the low pass filter implementation in the second Wikipedia article before with fairly good results, but I wasn't peak picking, just smoothing. I actually have a simple RC filter to electronically filter the data a little, but my electronics knowledge isn't quite up to understanding it properly. I'm not 100% sure that it's really helping at all, but at least it doesn't use up processor cycles (although, maybe it does lag the signal some, like I said, I don't understand it well enough).
|
|
|
|
|
Make several simple calculations and see which one works for you on sample/real data.
*slope / gradient (as mentioned in other replies)
*keep 6/10/12 etc.. Most recent highs - to get a rolling average of them.
*Exponential moving average (more responsive that a moving avergage)
Collect a few of these and you will find that one of them in comparison to the others gives you the indication than the samples are 'out of band'. Which is what I understand you to be seeking.
Most of the calculations above can be maintained on a rolling basis without having to keep too much data.
|
|
|
|
|
Try a running average of the last X samples and subtract it from the previous average. When the sign changes positive or negative, it's a peak or trough. You can change the value of X to filter out the noise from the patient moving.
|
|
|
|
|
General 
News 
Suggestion 
Question 
Bug 
Answer 
Joke 
Praise 
Rant 
Admin 
Submit an article or tip
