From fbce83293d5881cef3bbd0ce7ca3716dd8f0aad0 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Fri, 14 Nov 2014 07:40:00 +0100 Subject: [PATCH] =?UTF-8?q?initial=20import=20of=20classes=20for=20automat?= =?UTF-8?q?ic=20picking=20purposes=20[just=20imported=20by=20me;=20module?= =?UTF-8?q?=20has=20originally=20been=20written=20by=20Ludger=20K=C3=BCper?= =?UTF-8?q?koch]?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- pylot/core/pick/CharFuns.py | 342 ++++++++++++++++++++++++++++++++++++ 1 file changed, 342 insertions(+) create mode 100644 pylot/core/pick/CharFuns.py diff --git a/pylot/core/pick/CharFuns.py b/pylot/core/pick/CharFuns.py new file mode 100644 index 00000000..0fb5def1 --- /dev/null +++ b/pylot/core/pick/CharFuns.py @@ -0,0 +1,342 @@ +# -*- coding: utf-8 -*- +""" +Created Oct/Nov 2014 + +Implementation of the Characteristic Functions (CF) published and described in: + +Kueperkoch, L., Meier, T., Lee, J., Friederich, W., & EGELADOS Working Group, 2010: +Automated determination of P-phase arrival times at regional and local distances +using higher order statistics, Geophys. J. Int., 181, 1159-1170 + +Kueperkoch, L., Meier, T., Bruestle, A., Lee, J., Friederich, W., & EGELADOS +Working Group, 2012: Automated determination of S-phase arrival times using +autoregressive prediction: application ot local and regional distances, Geophys. J. Int., +188, 687-702. + +:author: MAGS2 EP3 working group +""" +import numpy as np +from obspy.core import Stream + +class CharacteristicFunction(object): + ''' + SuperClass for different types of characteristic functions. + ''' + def __init__(self, data, cut, t2, order, t1=None, fnoise=0.001): + ''' + Initialize data type object with information from the original + Seismogram. + + :param: data + :type: `~obspy.core.stream.Stream` + + :param: cut + :type: tuple + + :param: t2 + :type: float + + :param: order + :type: int + + :param: t1 + :type: float (optional, only for AR) + + :param: fnoise + :type: float (optional, only for AR) + ''' + + assert isinstance(data, Stream), "%s is not a Stream object" % str(data) + + self.orig_data = data[0] + self.dt = self.orig_data.stats.delta + self.setCut(cut) + self.setTime1(t1) + self.setTime2(t2) + self.setOrder(order) + self.setFnoise(fnoise) + self.calcCF(self.getDataArray()) + self.arpara = np.array([]) + self.xpred = np.array([]) + + def __str__(self): + return '''\n\t{name} object:\n + Cut:\t\t{cut}\n + t1:\t{t1}\n + t2:\t{t2}\n + Order:\t\t{order}\n + Fnoise:\t{fnoise}\n + '''.format(name=type(self).__name__, + cut=self.getCut(), + t1=self.getTime1(), + t2=self.getTime2(), + order=self.getOrder(), + fnoise=self.getFnoise()) + + def getCut(self): + return self.cut + + def setCut(self, cut): + self.cut = cut + + def getTime1(self): + return self.t1 + + def setTime1(self, t1): + self.t1 = t1 + + def getTime2(self): + return self.t2 + + def setTime2(self, t2): + self.t2 = t2 + + def getOrder(self): + return self.order + + def setOrder(self, order): + self.order = order + + def getIncrement(self): + return self.dt + + def getFnoise(self): + return self.fnoise + + def setFnoise(self, fnoise): + self.fnoise = fnoise + + def getCF(self): + return self.cf + + def getDataArray(self, cut=None): + ''' + If cut times are given, time series is cut from cut[0] (start time) + till cut[1] (stop time) in order to calculate CF for certain part + only where you expect the signal! + input: cut (tuple) () + cutting window + ''' + if cut is not None: + if self.cut[0] == 0: + start = 0 + else: + start = self.cut[0] / self.dt + stop = self.cut[1] / self.dt + data = self.orig_data.data[start:stop] + return data + + return self.orig_data.data + + def calcCF(self, data=None): + self.cf = data + + def arDet(self, data, order, rind, ldet): + pass + + def arPred(self, data, arpara, rind, lpred): + pass + + + +class AICcf(CharacteristicFunction): + ''' + Function to calculate the Akaike Information Criterion (AIC) after + Maeda (1985). + :param: data, time series (whether seismogram or CF) + :type: tuple + + Output: AIC function + ''' + + def calcCF(self, data): + print 'Calculating AIC ...' + xnp = self.getDataArray() + datlen = len(xnp) + k = np.arange(1, datlen) + cf = np.zeros(datlen) + cumsumcf = np.cumsum(np.power(xnp, 2)) + i = np.where(cumsumcf == 0) + cumsumcf[i] = np.finfo(np.float64).eps + cf[k] = ((k - 1) * np.log(cumsumcf[k] / k) + (datlen - k + 1) * + np.log((cumsumcf[datlen - 1] - + cumsumcf[k - 1]) / (datlen - k + 1))) + cf[0] = cf[1] + inf = np.isinf(cf) + ff = np.where(inf == 'True') + if len(ff) >= 1: + cf[ff] = 0 + + self.cf = cf - np.mean(cf) + + +class HOScf(CharacteristicFunction): + ''' + Function to calculate skewness (statistics of order 3) or kurtosis + (statistics of order 4), using one long moving window, as published + in Kueperkoch et al. (2010). + ''' + + def calcCF(self, data): + + xnp = self.getDataArray(self.getCut()) + + if self.getOrder() == 3: # this is skewness + print 'Calculating skewness ...' + y = np.power(xnp, 3) + y1 = np.power(xnp, 2) + elif self.getOrder() == 4: # this is kurtosis + print 'Calculating kurtosis ...' + y = np.power(xnp, 4) + y1 = np.power(xnp, 2) + + # Initialisation + # t2: long term moving window + ilta = round(self.getTime2() / self.getIncrement()) + lta = y[0] + lta1 = y1[0] + + # moving windows + LTA = np.zeros(len(xnp)) + for j in range(3, len(xnp)): + if j <= ilta: + lta = (y[j] + lta * (j - 1)) / j + lta1 = (y1[j] + lta1 * (j - 1)) / j + else: + lta = (y[j] - y[j - ilta]) / ilta + lta + lta1 = (y1[j] - y1[j - ilta]) / ilta + lta1 + # define LTA + if self.getOrder() == 3: + LTA[j] = lta / np.power(lta1, 1.5) + elif self.getOrder() == 4: + LTA[j] = lta / np.power(lta1, 2) + + LTA[0:3] = 0 + + self.cf = LTA + + +class ARZcf(CharacteristicFunction): + + def calcCF(self, data): + + print 'Calculating AR-prediction error from single trace ...' + xnp = self.getDataArray(self.getCut()) + + # some parameters needed + # add noise to time series + xnoise = xnp + np.random.normal(0.0, 1.0, len(xnp)) * self.getFnoise() * max(abs(xnp)) + tend = len(xnp) + # Time1: length of AR-determination window [sec] + # Time2: length of AR-prediction window [sec] + ldet = int(round(self.getTime1() / self.getIncrement())) # length of AR-determination window [samples] + lpred = int(np.ceil(self.getTime2() / self.getIncrement())) # length of AR-prediction window [samples] + + cf = [] + step = ldet + self.getOrder() - 1 + for i in range(ldet + self.getOrder() - 1, tend - lpred + 1): + if i == step: + ''' + In order to speed up the algorithm AR parameters are kept for time + intervals of length lpred + ''' + # determination of AR coefficients + self.arDet(xnoise, self.getOrder(), i - ldet, i) + step = step + lpred + + # AR prediction of waveform using calculated AR coefficients + self.arPred(xnp, self.arpara, i + 1, lpred) + # prediction error = CF + err = np.sqrt(np.sum(np.power(self.xpred[i:i + lpred] - xnp[i:i + lpred], 2)) / lpred) + cf.append(err) + + # convert list to numpy array + cf = np.asarray(cf) + self.cf = cf + + def arDet(self, data, order, rind, ldet): + ''' + Function to calculate AR parameters arpara after Thomas Meier (CAU), published + in Kueperkoch et al. (2012). This function solves SLE using the Moore- + Penrose inverse, i.e. the least-squares approach. + :param: data, time series to calculate AR parameters from + :type: array + + :param: order, order of AR process + :type: int + + :param: rind, first running summation index + :type: int + + :param: ldet, length of AR-determination window (=end of summation index) + :type: int + + Output: AR parameters arpara + ''' + + # recursive calculation of data vector (right part of eq. 6.5 in Kueperkoch et al. (2012) + rhs = np.zeros(self.getOrder()) + for k in range(0, self.getOrder()): + for i in range(rind, ldet): + rhs[k] = rhs[k] + data[i] * data[i - k] + + # recursive calculation of data array (second sum at left part of eq. 6.5 in Kueperkoch et al. 2012) + A = np.array([[0, 0], [0, 0]]) + for k in range(1, self.getOrder() + 1): + for j in range(1, k + 1): + for i in range(rind, ldet): + ki = k - 1 + ji = j - 1 + A[ki, ji] = A[ki, ji] + data[i - ji] * data[i - ki] + + A[ji, ki] = A[ki, ji] + + # apply Moore-Penrose inverse for SVD yielding the AR-parameters + self.arpara = np.dot(np.linalg.pinv(A), rhs) + + def arPred(self, data, arpara, rind, lpred): + ''' + Function to predict waveform, assuming an autoregressive process of order + p (=size(arpara)), with AR parameters arpara calculated in arDet. After + Thomas Meier (CAU), published in Kueperkoch et al. (2012). + :param: data, time series to be predicted + :type: array + + :param: arpara, AR parameters + :type: float + + :param: rind, first running summation index + :type: int + + :param: lpred, length of prediction window (=end of summation index) + :type: int + + Output: predicted waveform z + ''' + # be sure of the summation indeces + if rind < len(arpara) + 1: + rind = len(arpara) + 1 + if rind > len(data) - lpred + 1: + rind = len(data) - lpred + 1 + if lpred < 1: + lpred = 1 + if lpred > len(data) - 1: + lpred = len(data) - 1 + + z = np.append(data[0:rind], np.zeros(lpred)) + for i in range(rind, rind + lpred): + for j in range(1, len(arpara) + 1): + ji = j - 1 + z[i] = z[i] + arpara[ji] * z[i - ji] + + self.xpred = z + + +class ARHcf(CharacteristicFunction): + + pass + + +class AR3Ccf(CharacteristicFunction): + + pass