I goofed; peakdet is actually used

Author: pavelkomarov

pynumdiff/tests/test_utils.py

@@ -47,6 +47,34 @@ def test_convolutional_smoother():
     assert np.allclose(utility.convolutional_smoother(x, kernel_even, num_iterations=3), np.ones(len(x)))
 
 
+def test_peakdet(request):
+    """Verify peakdet finds peaks and valleys"""
+    t = np.arange(0, 10, 0.001)
+    x = 0.3*np.sin(t) + np.sin(1.3*t) + 0.9*np.sin(4.2*t) + 0.02*np.random.randn(10000)
+    maxtab, mintab = utility.peakdet(x, 0.5, t)
+
+    if request.config.getoption("--plot"):
+        from matplotlib import pyplot
+        pyplot.plot(t, x)
+        pyplot.plot(mintab[:,0], mintab[:,1], 'g*')
+        pyplot.plot(maxtab[:,0], maxtab[:,1], 'r*')
+        pyplot.title('peakdet validataion')
+        pyplot.show()
+
+    assert np.allclose(maxtab, [[0.447, 1.58575613], # these numbers validated by eye with --plot
+                                [1.818, 1.91349239],
+                                [3.316,-0.02740252],
+                                [4.976, 0.74512778],
+                                [6.338, 1.89861691],
+                                [7.765, 0.57577842],
+                                [9.402, 0.59450898]])
+    assert np.allclose(mintab, [[1.139, 0.31325728],
+                                [2.752,-1.12769567],
+                                [4.098,-2.00326946],
+                                [5.507,-0.31714122],
+                                [7.211,-0.59708324],
+                                [8.612,-1.7118216 ]])
+
 def test_slide_function():
     """Verify the slide function's weighting scheme calculates as expected"""
     def identity(x, dt): return x, 0 # should come back the same

pynumdiff/utils/utility.py

@@ -130,3 +130,54 @@ def slide_function(func, x, dt, kernel, *args, stride=1, pass_weights=False, **k
         weight_sum[window] += w # save sum of weights for normalization at the end
 
     return x_hat/weight_sum, dxdt_hat/weight_sum
+
+
+def peakdet(x, delta, t=None):
+    """Find peaks and valleys of 1D array. A point is considered a maximum peak if it has the maximal
+    value, and was preceded (to the left) by a value lower by delta. Converted from MATLAB script at
+    http://billauer.co.il/peakdet.html Eli Billauer, 3.4.05 (Explicitly not copyrighted). This function
+    is released to the public domain; Any use is allowed.
+
+    :param np.array[float] x: array for which to find peaks and valleys
+    :param float delta: threshold for finding peaks and valleys. A point is considered a maximum peak
+        if it has the maximal value, and was preceded (to the left) by a value lower by delta.
+    :param np.array[float] t: optional domain points where data comes from, to make indices into locations
+
+    :return: tuple[np.array, np.array] of\n
+             - **maxtab** -- indices or locations (column 1) and values (column 2) of maxima
+             - **mintab** -- indices or locations (column 1) and values (column 2) of minima
+    """
+    maxtab = []
+    mintab = []
+    if t is None:
+        t = np.arange(len(x))
+    elif len(x) != len(t):
+        raise ValueError('Input vectors x and t must have same length')
+    if not (np.isscalar(delta) and delta > 0):
+        raise ValueError('Input argument delta must be a positive scalar')
+
+    mn, mx = np.inf, -1*np.inf
+    mnpos, mxpos = np.nan, np.nan
+    lookformax = True
+    for i in np.arange(len(x)):
+        this = x[i]
+        if this > mx:
+            mx = this
+            mxpos = t[i]
+        if this < mn:
+            mn = this
+            mnpos = t[i]
+        if lookformax:
+            if this < mx-delta:
+                maxtab.append((mxpos, mx))
+                mn = this
+                mnpos = t[i]
+                lookformax = False # now searching for a min
+        else:
+            if this > mn+delta:
+                mintab.append((mnpos, mn))
+                mx = this
+                mxpos = t[i]
+                lookformax = True # now searching for a max
+
+    return np.array(maxtab), np.array(mintab)