From bb84e27e54b3bc54d31809a4f7a2fa96c19e6b75 Mon Sep 17 00:00:00 2001
From: Sebastian Wehling-Benatelli <sebastian.wehling@rub.de>
Date: Tue, 21 Apr 2015 08:14:26 +0200
Subject: [PATCH] implemented a much faster calculation of the zero crossing
 beside an average calculation over the whole signal window of the dominant
 period

---
 pylot/core/pick/utils.py | 38 ++++++++++++++++----------------------
 1 file changed, 16 insertions(+), 22 deletions(-)

diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py
index e97dab3f..9cdacb49 100644
--- a/pylot/core/pick/utils.py
+++ b/pylot/core/pick/utils.py
@@ -53,32 +53,22 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
     #get signal window
     isignal = getsignalwin(t, Pick1, TSNR[2])
     #calculate noise level
-    nlevel = max(abs(x[inoise])) * nfac
+    nlevel = np.sqrt(np.mean(np.square(x[inoise]))) * nfac
     #get time where signal exceeds nlevel
-    ilup = np.where(x[isignal] > nlevel)
-    ildown = np.where(x[isignal] < -nlevel)
-    if len(ilup[0]) <= 1 and len(ildown[0]) <= 1:
-        print 'earllatepicker: Signal lower than noise level, misspick?'
-        return
-    il = min([ilup[0][0], ildown[0][0]])
+    ilup, = np.where(x[isignal] > nlevel)
+    ildown, = np.where(x[isignal] < -nlevel)
+    if not ilup.size and not ildown.size:
+        raise ValueError('earllatepicker: Signal lower than noise level')
+    il = min(np.min(ilup) if ilup.size else float('inf'),
+             np.min(ildown) if ildown.size else float('inf'))
     LPick = t[isignal][il]
 
     #get earliest possible pick
-    #get next 2 zero crossings after most likely pick
-    #initial onset is assumed to be the first zero crossing
-    zc = []
-    zc.append(Pick1)
-    i = 0
-    for j in range(isignal[0][1], isignal[0][len(t[isignal]) - 1]):
-        i = i + 1
-        if x[j - 1] <= 0 and x[j] >= 0:
-            zc.append(t[isignal][i])
-        elif x[j - 1] > 0 and x[j] <= 0:
-            zc.append(t[isignal][i])
-        if len(zc) == 3:
-            break
-    #calculate maximum period T0 of signal out of zero crossings
-    T0 = max(np.diff(zc))  #this is half wave length!
+
+    #determine all zero crossings in signal window
+    zc = crossings_nonzero_all(x[isignal])
+    #calculate mean half period T0 of signal as the average of the
+    T0 = np.mean(np.diff(zc)) * X[0].stats.delta  #this is half wave length!
     #T0/4 is assumed as time difference between most likely and earliest possible pick!
     EPick = Pick1 - T0 / 2
 
@@ -288,6 +278,10 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
 
     return FM
 
+def crossings_nonzero_all(data):
+    pos = data > 0
+    npos = ~pos
+    return ((pos[:-1] & npos[1:]) | (npos[:-1] & pos[1:])).nonzero()[0]
 
 def getSNR(X, TSNR, t1):
     '''