sandell@ferret.ils.nwu.edu (Greg Sandell) (10/24/90)
------------------------------------------------------------------
A few weeks back I posted a query about peak detection algorithms
and have since recieved a number of replies. A summary of the
replies appears below, following the original query I posted.
------------------------------------------------------------------
[ Original posting begins here ]
Can anybody tell me about the existence of any peak detection
algorithms? Below is an example of what I mean:
B & E are real peaks, the rest are not.
|--------------------------------------------------|
| B |
| * C |
| * * * |
| * ** * E |
| * * * * |
| * * * ** |
| * * * * |
| * * * * |
| * * D * * |
| * * * * * |
| * * ** * * |
| A * * * * * |
|-------*-----------------*-----*--------------*---| (threshold)
| * * * * * |
| * * * * * * |
|* * * * *|
|* * *|
|--------------------------------------------------|
A is not a real peak because it falls below my arbitrary
threshold. C is not real because it is subsidiary to
the fall away from real peak B, and D is not real
because it is subsidiary to the climb toward real
peak E.
Certainly any number of hacks will solve this (with
variables such as threshold, minimum valley required
before and after a peak, and so on). But I would
like to know if there are any classic solutions to this
problem. Also, has code for such an algorithm been
printed anywhere (such as Numerical Recipes)?
To satsify your curiosity, I need a peak detector for
a speech processing problem: I have a curve which
represents the gross changes in amplitude for a recording
of speech. By identifying the primary peaks, I can
find the principle syllable boundaries. Obviously this
type of signal proposes particular issues for peak
detection: for example, the peaks are skewed to the right
in speech (rising side is more acute than the falling side),
because onsets are quicker than decays.
I would appreciate replies via email. I will be glad to share
any code I receive with all who are interested.
Thanks in advance,
Greg Sandell
****************************************************************
* Greg Sandell (sandell@ils.nwu.edu) Evanston, IL *
* Institute for the Learning Sciences, Northwestern University *
****************************************************************
[ End of original posting ]
------------------------------------------------------------------
A few people told me that they were also working on problems that
required peak detection algorithms. The variety of problems they
described underscored the fact that every application involves
its own idiosyncracies, as indicated below:
------------------------------------------------------------------
Julian Vrieslander eacj@theory.tn.cornell.edu wrote:
> I need to be able to
> determine the time-of-occurrence of neural spikes in a stream of samples
> coming from an microelectrode->amplifier->ADC circuit.
Robert John Butera Jr. gt5614b@prism.gatech.edu wrote:
> In a class
> I took recently which spent three weeks in radar processing, we had to pick
> the N largest peaks from a signal (similar to what you are trying to do)
> and avoid any "sub" or "transitional" peaks.
sinistar@mrcnext.cso.uiuc.edu (Jeff O'Hare) wrote:
> This application is for checking major flooding events...it uses certain
> criterion that only could apply to river flow data.
rml@extro.ucc.su.OZ.AU (Dr.Ross Lazarus) wrote:
> there's a whole raftload of work in this area in endocrine physiology -
> most endocrine signalling is done with pulsatile patterns of release.
> The problem is complicated by the high noise levels due to the analytical
> methods for hormone detection - ria.
------------------------------------------------------------------
There were, however, a few general solutions that arose from this query.
First let us define simplepeaks(), which is a poor peak-detector, but
which will come in handy for post-processing. This simply looks for any
triad of data points which form a "peak-like" or ^ shape where the middle
point is above threshold. (Apologies in advance to people who don't like
C code!)
------------------------------------------------------------------
int simplepeaks(data,N,threshold,peaktimes)
int N, *peaktimes;
double *data, threshold;
{
int i, peakcount;
double left, middle, right;
for (i = 1, peakcount = 0; i < N-2; i++) {
left = data[i-1]; /* left, middle and right used here */
middle = data[i]; /* to make strategy clear. */
right = data[i+1];
if ((middle > threshold) && (left < middle) && (right < middle)) {
peaktimes[peakcount++] = i;
}
}
return(peakcount);
}
------------------------------------------------------------------
The peak detector which worked best for my problem was a simple lowpass
filter, combined with simplepeaks(). A lowpass filter treats smaller,
subsidiary bumps as high frequency information to be filtered out, so
(returning to my original example) B & C get smoothed down to one peak,
as does D & E. My application involved the analysis of spoken speech,
which tends to have a rate of 200-300 syllables a minutes, that is,
3-5 Hz. So I could use the lowpass filter to get rid of anything above
5 Hz or so. Even if another application was not so periodic as speech,
it would always be possible to run an FFT on the data first to find
out what range of frequencies were involved, then make a guess as to
what frequencies to cut off. Some code for this is shown below.
"srate" means sampling rate, the number of values per second. If
an application is not time-dependent, you can select an arbitrary
value, but then you have to find a meaningful cutoff_freq.
------------------------------------------------------------------
#define PI 3.141592564
#include <math.h>
void bpass_array(data,N,srate,cutoff_freq)
double *data, srate, cutoff_freq;
int N;
{
int i;
double high, low, range, fsample;
double cfreq, in, theta, omega, r, b1, b2, k;
double old1, old2; /* these are all bandpass filter related */
cfreq = cutoff_freq/2.0;
theta = 2.0 * PI * (cfreq/srate);
omega = 2.0 * PI * (cutoff_freq/srate);
r = sqrt(exp(-1.0 * omega));
b1 = 2.0 * r * cos(theta);
b2 = -1.0 * (r * r);
k = (1.0-r)*sqrt(1.0+(r * r)+(2.0*r)-(4.0*r*(cos(theta) * cos(theta))));
old1 = old2 = 0.0;
/* My data consisted of 16 bit audio samples, meaning they varied
from -32767 to +32767. My bandpass filter, however, wants to
see this scaled down to values between -1.0 and +1.0 */
for (i = 0, low = 40000.0, high = -40000.0; i < N; i++) {
high = data[i] > high ? data[i] : high;
low = data[i] < low ? data[i] : low;
}
range = high - low;
for (i = 0; i < N; i++) {
/* scales data down to -1 to +1, does the lowpass */
fsample = (((data[i]-low)/range) * k) + (old1 * b1) + (old2 * b2);
old2 = old1;
old1 = fsample;
/* scales it back to the original range */
data[i] = (fsample * range) + low;
}
}
------------------------------------------------------------------
The data that I included in my original posting turns into this when
I run it through bpass_array():
------------------------------------------------------------------
|--------------------------------------------------|
| 1112 |
| 11 1 |
| 1 1 |
| 1 1 |
| 1 |
| 1 1 2 11 |
| 1 1 1 |
| 1 1 1 |
| 1 1 1 |
| 1 |
| 1 1 1 1|
| 1 1 1|
| 1 1 11 |
| 1 11 |
| 1 1111 |
| 2 |
| 11 |
| 11 |
|211 |
|--------------------------------------------------|
1.000 58.000
------------------------------------------------------------------
This transforms the data to something that makes the peaks easy
to find, such as by the simplepeaks() function.
The trouble with bpass_array() is that the location of the peak has been
generalized to a particular region; although this was fine for my
application, some people will find this too imprecise (i.e. they
want to know exactly which datapoint in the *original* data is the peak).
In that case, an extra function will be needed to correct the error.
Robert John Butera Jr. gt5614b@prism.gatech.edu applied the lowpass
filter to his radar data (see above), and here's his report:
------------------------------------------------------------------
> The class agreed that there
> was no known "elegant" way to do this, at least to us. The method I used
> required some knowledge of the characteristics of the input data.
> A "heavy" but quick (second order) low pass filter was applied to the data.
> The highest point of the sample spectrum was selected as the first peak.
> Then the area around the peak (an arbtrarily set number of samples) were
> reset to values below the threshold, and the process repeated for the next peak.
>
> This is obviously a "klugy" process, since we had a finite set of data and
> knew about how many peaks to look for. Good luck!
------------------------------------------------------------------
I told Carsten Mandal <cfm@imada.dk> about my lowpass solution and he
(she?) sent the following comment:
------------------------------------------------------------------
> Yes, that is true - but when it comes to data, where you have
> a higher rate, that you wish to detect i.e. that there is not much
> difference between "the good vs. the bad" freqs. you are into trouble -
> at least that's what I have experienced. When it comes to detecting
> a low volume fixed pitch tone mixed with higher leveled speek, this
> method have some drawbacks - but making a bandpass (3. order) filter
> and making some statictics on the output from the filter seems to
> solve that problem.
>
> Thanks for the answer.
>
> Carsten Mandal
------------------------------------------------------------------
black@seismo.CSS.GOV (Mike Black) sent some elegant code
in C which defined a peak as the highest above-threshold point laying
between two drops to below-threshold. This solution was not so
great for my data because there were multiple peaks inbetween drops
to threshold that I wanted to save. However, this approach is great
for cleaning up data that looks like this:
------------------------------------------------------------------
|--------------------------------------------------|
| 1 |
| 1 11 1 |
| |
| 11 1 1 111 |
| 1 1 111 1 1 1 1 1 |
| |
| 1 11 1 1 111 111 |
| 1 1 1 1 1 1 1 |
| |
| 1 1 1 1 |
| |
| 1 1 1 1 |
| 1 1 1 1 |
| |
| 1 1 1 1 |
| 1 1 1 1 |
| |
| 1 1 1 1 |
|22 1211 12122|
|--------------------------------------------------|
1.000 81.000
------------------------------------------------------------------
The code will return the values and locations of the topmost
values of these two peaks.
Here is Mike Black's code, with some minor improvements from me:
------------------------------------------------------------------
peaks(x, N, THRESHOLD)
int *x, N, THRESHOLD;
{
int i,above_threshold=0,peak_loc=-1;
int max=0;
for (i=0;i<N-1;i++) {
if (x[i] >= THRESHOLD) {
above_threshold = 1;
if (x[i] > max) {
max = x[i];
peak_loc = i;
}
}
else {
if (above_threshold) {
printf("Peak is at %d = %d\n",peak_loc,max);
above_threshold = 0;
}
peak_loc=-1;
max=0;
}
if (above_threshold) {
printf("Possible peak at %d\n",peak_loc);
}
}
------------------------------------------------------------------
After I told Mike Black about my success in using a lowpass filter, he
sent me the following alternative solution. I include it without any
further comment since I never investigated this approach.
------------------------------------------------------------------
From black@seismo.CSS.GOV (Mike Black):
> You could also use a weighted differential like this:
> (This assumes your peaks have some width which my other suggestion did not)
>
> Create a filter kernel like this:
>
> int kernel[] = { -1,-1,-1,-1,0,1,1,1,1 };
>
> Notice you have negative weights on one side and positive on the other.
> You can adjust the width of this to your data (more pairs of -1,1).
>
> The effect this has is just like a lowpass filter, but it also removes the
> DC offset. You apply this kernel at each point in your waveform.
>
> for (j=0;i<N-sizeof(kernel);j++)
> for (i=0;i<9,i++) result[j] += wave[j]*kernel[i];
> /* wave should be normalized (divide by 8) at this point, but it's unnecessary
>
> This can also be expanded out for speed:
>
> result[j] = - wave[j] - wave[j+1] - wave[j+2] - wave[j+3] + wave[j+5]
> + wave[j+6] + wave[j+7] + wave[j+8]
>
> You lose a few points at the beginning of the waveform doing this, but this
> is fairly fast. Your peaks show up at zero crossing points. All you need to
> do is apply a threshold against the result and you can find the peaks real
> easy.
------------------------------------------------------------------
rml@extro.ucc.su.OZ.AU (Dr.Ross Lazarus), who was dealing with endocrine
pulses (see above), sent the following elaborative information:
------------------------------------------------------------------
> The 2 most useful starting points would be
> 1) Veldhuis and Johnson, "Cluster analysis : a simple versatile and robust
> algorithm for endocrine pulse detection", Am.J.Physiol., 250 (Endocrinol.
> Metab 13):E486-493,1986
> 2) Merriam and Wachter,"Algorithms for the study of episodic hormone secretion",
>
> Am.J.Physiol. 243 (Endocrinol. Metab 6):E310-E318,1982
>
> There are no assumption free methods (!) yet widely known. They all make some
> clumsy requirement for the user to specify some arbitrary parameters. Veldhuis'
> method requires "t" statistics and cluster sizes. It suffers very badly from
> edge effects. Merriam's method requires specifying "G" parameters which are
> in practice something like t statistics.
>
> Veldhuis takes groups of points (clusters - of a size the user
> specifies for peaks and valleys or nadirs as he calls them) and compares their
> values with the preceding group using a t statistic - but there's a significant
> wrinkle - in the measurement of endocrine hormones at extremely small levels in
> blood, radio-immune assay is used and each assay has a measurable precision
> profile which is often formulated as a relationship between the sd and
> the mean of replicates. As you can see, it's starting to get complex ! The
> relationship between mean value and sd is used in this t testing to detect
> "significant" rises and falls. A peak is then defined as a significant rise
> preceeded and followed by a significant fall - hence the edge effect
> problem.
> Merriam's method is (in my view) nicer. He uses LOWESS smoothing to estimate
> a baseline through the data, then looks for significant rises with different
> (usually falling) levels regarded as significant with increasing numbers of
> points under examination - the window is normally limited to 5 points. When a
> significant point or group of points is found, these are downweighted in
> another iteration of the lowess smoothing - so the baseline goes down and
> another cycle of peak detection follows. Iterate until no more peaks are
> found and you've got PULSAR. Again, knowledge of the assay precision
> profile is used to estimate the size of the rises in terms of z scores.
> Your problem as I recall concerns sounds. Measurement is nothing like as
> problematic as ria where a 15% cv is considered precise !
------------------------------------------------------------------------
Well, that's it. Thanks again to the people who sent replies, esp.
to Mike Black for his helpful code.
Greg Sandell
------------------------------------------------------------------------
****************************************************************
* Greg Sandell (sandell@ils.nwu.edu) Evanston, IL *
* Institute for the Learning Sciences, Northwestern University *
****************************************************************