|
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.
|
|
|
|
|
I've done this with heart rate variability (HRV) and it's rather simple. Determine the frequency range of your useful data and setup a bandpass IIR or FIR filter to filter out all unwanted noise. Use matlab or similar tool to determine the coefficients of the filter. IIR filters are very easy to implement in code as they are just some floating point multiply ops. After your data is filtered, peaks can be easily identified with a simple rule e.g. if(sample[t] < sample[t-1]). All is in the design of the filter (compromise between delay vs and bandwidth).
modified on Thursday, June 2, 2011 3:38 AM
|
|
|
|
|
I found the following article useful when I was trying to find peaks & troughs in some data:-
http://billauer.co.il/peakdet.html[^]
Looking at your graph the only slight hiccup would be at 31s, but you can set the delta value to an appropriate value to overcome this I think.
|
|
|
|
|
Hi,
One good way to remove spikes is to use a median filter. Keep a window of the data, sort the values in the window, and use the one in the center of the range.
Don't know if that is appropriate for a real time application, but we use it for seismic data post-processing all the time.
Is that sleep-study data? I had to have one a few weeks ago, I have sleep apnea!
Regards,
Tom
|
|
|
|
|
I'm sure there's a relatively-straightforward way to do this, but my research so far has been a bit overwhelming... Heavy statistics and calculus have never been one of my strengths...
Basically, I need to calculate a best-fit curve on a scatter graph... Points are not evenly spaced on either axis, and there's no cost function to work with (These are derivations of stock/bond price histories)... Just raw (X,Y) values.
I've been trying to work out an algorithm that'll start with a second-degree polynomial function, iterate through and minimize the difference, and then jump to third-degree if it's not within tolerance (And so on, up to a maximum degree). Juggling the data around is something I can do blindfolded, but there's one thing I can't figure out...
Given the polynomial coefficients ((a,b,c,d) -> y = a + bx + cx^2 + dx^3) and a function to determine the total difference across the curve (Or the total square difference, or whatever's needed), how do I adjust the coefficients for the next pass? I mean, how do I choose new values?
I was originally going to try to fix the highest-degree coefficient first, then drill down to the zeroth-degree one (I have a root-solving algorithm that I use for yield calculations - Was going to adapt that), but then I looked at it again with the proper amount of caffeine, and realized that solving one coefficient at a time wasn't going to work
SOLVED: See my reply below... Linear/Polynomial Regression
|
|
|
|
|
Not sure if this will help but, have you tried a spline interpolation algorithm?
|
|
|
|
|
Kind of the opposite of where I'm going here... The goal isn't to pass through every point, but to give more of a trend-line.
Generally going to be working with a few hundred data points, but only looking for maybe a second or third-degree function in most cases.
<small><i>Proud to have finally moved to the A-Ark. <a href="http://en.wikipedia.org/wiki/Places_in_The_Hitchhiker%27s_Guide_to_the_Galaxy#Golgafrincham">Which one are you in?</a></i><br>Author of the <a href="http://www.serpentooth.com">Guardians Saga</a> (Sci-Fi/Fantasy novels)</br></small>
|
|
|
|
|
Ok, now I understand it better. Then maybe bezier curves would be a better choice.
|
|
|
|
|
Hi Ian,
I don't think there are any closed formula's for what you are trying. What I would recommend is this, inspired by Newton-Raphson:
1. have some values a1,b1,c1,d1 and the corresponding value of your error function e1.
2. slightly change a1 (say by 1%) and recalculate e; find the "a-sensitivity" as (change in e)/(change in a)
3. same for b1, then c1, then d1
4. now choose new values a2,b2,c2,d2 based on the four sensitivities, maybe like so: a2=a1+coef*e/"a-sensitivity" where maybe coef is one.
then iterate until convergence is reached (i.e. sufficiently accurate); warning: this might diverge right away, if so try smaller values of coef (you could halve it in each consecutive try).
FWIW: with more-or-less random data, you'll never find a good fit with a single polynomial. What one normally does is a local approximation, which applies in a small region only. That is similar to a font being drawn as a sequence of Bezier curves; it pays of to have several simple Bezier curves, rather than trying to have a single very complex one (which in general will not be good enough).
Luc Pattyn [Forum Guidelines] [My Articles] Nil Volentibus Arduum
Please use <PRE> tags for code snippets, they preserve indentation, improve readability, and make me actually look at the code.
|
|
|
|
|
Hmm, didn't think of the sensitivity thing... Giving that a shot now.
FWIW: with more-or-less random data, you'll never find a good fit with a single polynomial. What one normally does is a local approximation, which applies in a small region only. That is similar to a font being drawn as a sequence of Bezier curves; it pays of to have several simple Bezier curves, rather than trying to have a single very complex one (which in general will not be good enough).
It's not actually random data... It's graphing related values, so it should roughly cluster around y=x unless there's another factor (Which is what I'm trying to illustrate).
Proud to have finally moved to the A-Ark. Which one are you in?
Author of the Guardians Saga (Sci-Fi/Fantasy novels)
|
|
|
|
|
General 
News 
Suggestion 
Question 
Bug 
Answer 
Joke 
Praise 
Rant 
Admin 
Submit an article or tip
