Improving wave detection algorithm

Question

I am new to python and programming and working on the wave peak detection algorithm on my spline - interpolated plot. I have used the code given on this link : https://gist.github.com/endolith/250860 . I have to make this algorithm work on any type of waveform i.e. low as well as high amplitude, baseline not aligned, etc. The goal is to calculate the no of waves in the plot. But my peak detection calculates "invalid" peaks and so gives the wrong answers. By "invalid" peaks I mean if there are two notches close to each other at the wave peak, the program detects 2 peaks i.e. 2 waves when actually its just 1 wave. I have tried changing the 'delta' parameter defined in the peak detection function given on the link but that doesn't solve the generalisation goal which I am working on..Please suggest any improvement on the algorithm or any other approach which I should be using. Any kind of help is welcomed. Thanks in advance.

P.S. I am unable to upload an image of the wrong peak-detected wave plot. I hope my explanation is sufficient enough... Code is as follows

wave = f1(xnew)/(max(f1(xnew))) ##interpolated wave
maxtab1, mintab1 = peakdet(wave,.005) 
maxi = max(maxtab1[:,1])
for i in range(len(maxtab1)):
    if maxtab1[i,1] > (.55 * maxi) : ## Thresholding
        maxtab.append((maxtab1[i,0],maxtab1[i,1]))
arr_maxtab = np.array(maxtab)
dist = 1500 ## Threshold defined for minimum distance between two notches to be considered as two separate waves
mtdiff = diff(arr_maxtabrr[:,0])
final_maxtab = []
new_arr = []
for i in range(len(mtdiff)):
if mtdiff[i] < dist :
            new_arr.append((arr_maxtab[i+1,0],arr_maxtab[i+1,1])) 
for i in range(len(arr_maxtab)):
if  not arr_maxtab[i,0] in new_arr[:,0]:
    final_maxtab.append((arr_maxtab[i,0], arr_maxtab[i,1]))
else:
final_maxtab = maxtab

Answer 1

The ability to distinguish notches from true peaks implies you have a fundamental spacing between peaks. Said differently, there is a minimum frequency resolution at which you would like to run your peak detection search. If you zoom into a signal at which you are well narrower than the noise floor, you'll observe zig zags that seem to 'peak' every few samples.

What it sounds like you want to do is the following:

1. Smooth the signal.
2. Find the 'real' peaks.

Or more precisely,

1. Run the signal through a low pass filter.
2. Find the peaks within your acceptable peak widths with sufficient signal to noise ratio.

Step 1: Low-Pass Filtering

To do step 1, I recommend you use the signal processing tools provided by scipy.

I'm adapted this cookbook entry, which shows how to use FIR filters to lowpass a signal using scipy.

from numpy import cos, sin, pi, absolute, arange, arange
from scipy.signal import kaiserord, lfilter, firwin, freqz, find_peaks_cwt
from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, grid, axes, show, scatter

#------------------------------------------------
# Create a signal. Using wave for the signal name.
#------------------------------------------------

sample_rate = 100.0
nsamples = 400
t = arange(nsamples) / sample_rate
wave = cos(2*pi*0.5*t) + 0.2*sin(2*pi*2.5*t+0.1) + \
        0.2*sin(2*pi*15.3*t) + 0.1*sin(2*pi*16.7*t + 0.1) + \
            0.1*sin(2*pi*23.45*t+.8)


#------------------------------------------------
# Create a FIR filter and apply it to wave.
#------------------------------------------------

# The Nyquist rate of the signal.
nyq_rate = sample_rate / 2.0

# The desired width of the transition from pass to stop,
# relative to the Nyquist rate.  We'll design the filter
# with a 5 Hz transition width.
width = 5.0/nyq_rate

# The desired attenuation in the stop band, in dB.
ripple_db = 60.0

# Compute the order and Kaiser parameter for the FIR filter.
N, beta = kaiserord(ripple_db, width)

# The cutoff frequency of the filter.
cutoff_hz = 10.0

# Use firwin with a Kaiser window to create a lowpass FIR filter.
taps = firwin(N, cutoff_hz/nyq_rate, window=('kaiser', beta))

# Use lfilter to filter x with the FIR filter.
filtered_x = lfilter(taps, 1.0, wave)

Step 2: Peak Finding

For step 2, I recommend you use the wavelet transformation peak finder provided by scipy. You must provide as an input your filtered signal and a vector running from the minimum to maximum possible peak widths. This vector will be used as the basis of the wavelet transformation.

#------------------------------------------------
# Step 2: Find the peaks
#------------------------------------------------

figure(4)
plot(t[N-1:]-delay, filtered_x[N-1:], 'g', linewidth=1)

peakind = find_peaks_cwt(filtered_x, arange(3,20))
scatter([t[i] - delay for i in peakind], [filtered_x[i] for i in peakind], color="red")
for i in peakind:
    print t[i] + delay

xlabel('t')
grid(True)

show()