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