From 45eb138717a1cf4b6ece730ef32a669c32116cae Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= <kueperkoch@bestec-for-nature.com>
Date: Mon, 30 Mar 2015 10:58:24 +0200
Subject: [PATCH] New function to derive SNR, returns SNR and SNR[dB].

---
 pylot/core/pick/getSNR.py | 66 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 66 insertions(+)
 create mode 100644 pylot/core/pick/getSNR.py

diff --git a/pylot/core/pick/getSNR.py b/pylot/core/pick/getSNR.py
new file mode 100644
index 00000000..6c2b6872
--- /dev/null
+++ b/pylot/core/pick/getSNR.py
@@ -0,0 +1,66 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+   Created Mar/Apr 2015
+   Function to calculate SNR of certain part of seismogram relativ 
+   to given time. Returns SNR and SNR [dB].
+
+   :author: Ludger Kueperkoch /MAGS EP3 working group
+"""
+from obspy.core import Stream
+import numpy as np
+
+def getSNR(X, TSNR, t1):
+    '''
+    Function to calculate SNR of certain part of seismogram relative to
+    given time (onset) out of given noise and signal windows. A safety gap
+    between noise and signal part can be set. Returns SNR and SNR [dB].
+
+    :param: X, time series (seismogram)
+    :type:  `~obspy.core.stream.Stream`
+
+    :param: TSNR, length of time windows [s] around t1 (onset) used to determine SNR
+    :type: tuple (T_noise, T_gap, T_signal)
+
+    :param: t1, initial time (onset) from which noise and signal windows are calculated
+    :type: float
+    '''
+
+    assert isinstance(X, Stream), "%s is not a stream object" % str(X)
+
+    SNR = None
+    SNRdB = None
+    x = X[0].data
+    t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, X[0].stats.delta)
+    #some parameters needed:
+    tnoise = TSNR[0]  #noise window length for calculating noise level
+    tsignal = TSNR[2] #signal window length
+    tsafety = TSNR[1] #safety gap between signal onset and noise window
+    #get noise window
+    inoise = np.where((t <= max([t1 - tsafety, 0])) \
+             & (t >= max([t1 - tnoise - tsafety, 0])))
+    #get signal window
+    isignal = np.where((t <= min([t1 + tsignal + tsafety, len(x)])) \
+              & (t >= t1))
+    if np.size(inoise) < 1:
+       print 'getSNR: Empty array inoise, check noise window!'
+       return
+    elif np.size(isignal) < 1:
+       print 'getSNR: Empty array isignal, check signal window!'
+       return
+
+    #calculate ratios
+    SNR = max(abs(x[isignal])) / np.mean(abs(x[inoise]))
+    SNRdB = 20 * np.log10(SNR)
+    
+    return SNR, SNRdB
+
+if __name__ == "__main__":
+   parser = argparse.ArgumentParser()
+   parser.add_argument('--X', type=~obspy.core.stream.Stream, help='time series (seismogram) read with obspy module read')
+   parser.add_argument('--TSNR', type=tuple, help='length of time windows around pick used to determine SNR \
+                       [s] (Tnoise, Tgap, Tsignal)')
+   parser.add_argument('--t1', type=float, help='initial time from which noise and signal windows are calculated')
+   args = parser.parse_args()
+   getSNR(args.X, args.TSNR, args.t1)
+