1# -*- coding: utf-8 -*-
2"""
3Signal processing functions for RtMemory objects.
4
5For sequential packet processing that requires memory (which includes recursive
6filtering), each processing function (e.g., :mod:`obspy.realtime.signal`)
7needs to manage the initialization and update of
8:class:`~obspy.realtime.rtmemory.RtMemory` object(s), and needs to know when
9and how to get values from this memory.
10
11For example: Boxcar smoothing: For each new data point available past the end
12of the boxcar, the original, un-smoothed data point value at the beginning of
13the boxcar has to be subtracted from the running boxcar sum, this value may be
14in a previous packet, so has to be retrieved from memory see
15:func:`obspy.realtime.signal.boxcar`.
16
17:copyright:
18    The ObsPy Development Team (devs@obspy.org), Anthony Lomax & Alessia Maggi
19:license:
20    GNU Lesser General Public License, Version 3
21    (https://www.gnu.org/copyleft/lesser.html)
22"""
23from __future__ import (absolute_import, division, print_function,
24                        unicode_literals)
25from future.builtins import *  # NOQA
26
27import math
28import sys
29
30import numpy as np
31
32from obspy.core.trace import Trace, UTCDateTime
33from obspy.realtime.rtmemory import RtMemory
34
35
36_PI = math.pi
37_TWO_PI = 2.0 * math.pi
38_MIN_FLOAT_VAL = 1.0e-20
39
40
41def offset(trace, offset=0.0, rtmemory_list=None):  # @UnusedVariable
42    """
43    Add the specified offset to the data.
44
45    :type trace: :class:`~obspy.core.trace.Trace`
46    :param trace: :class:`~obspy.core.trace.Trace` object to append to this
47        RtTrace
48    :type offset: float, optional
49    :param offset: offset (default is 0.0)
50    :type rtmemory_list: list of :class:`~obspy.realtime.rtmemory.RtMemory`,
51        optional
52    :param rtmemory_list: Persistent memory used by this process for specified
53        trace
54    :rtype: NumPy :class:`numpy.ndarray`
55    :return: Processed trace data from appended Trace object
56    """
57
58    if not isinstance(trace, Trace):
59        msg = "Trace parameter must be an obspy.core.trace.Trace object."
60        raise ValueError(msg)
61
62    trace.data += offset
63    return trace.data
64
65
66def scale(trace, factor=1.0, rtmemory_list=None):  # @UnusedVariable
67    """
68    Scale array data samples by specified factor.
69
70    :type trace: :class:`~obspy.core.trace.Trace`
71    :param trace:  :class:`~obspy.core.trace.Trace` object to append to this
72        RtTrace
73    :type factor: float, optional
74    :param factor: Scale factor (default is 1.0).
75    :type rtmemory_list: list of :class:`~obspy.realtime.rtmemory.RtMemory`,
76        optional
77    :param rtmemory_list: Persistent memory used by this process for specified
78        trace.
79    :rtype: NumPy :class:`numpy.ndarray`
80    :return: Processed trace data from appended Trace object.
81    """
82    if not isinstance(trace, Trace):
83        msg = "trace parameter must be an obspy.core.trace.Trace object."
84        raise ValueError(msg)
85    # XXX not sure how this should be for realtime analysis, here
86    # I assume, we do not want to change the underlying dtype
87    trace.data *= np.array(factor, dtype=trace.data.dtype)
88    return trace.data
89
90
91def integrate(trace, rtmemory_list=None):
92    """
93    Apply simple rectangular integration to array data.
94
95    :type trace: :class:`~obspy.core.trace.Trace`
96    :param trace:  :class:`~obspy.core.trace.Trace` object to append to this
97        RtTrace
98    :type rtmemory_list: list of :class:`~obspy.realtime.rtmemory.RtMemory`,
99        optional
100    :param rtmemory_list: Persistent memory used by this process for specified
101        trace.
102    :rtype: NumPy :class:`numpy.ndarray`
103    :return: Processed trace data from appended Trace object.
104    """
105    if not isinstance(trace, Trace):
106        msg = "trace parameter must be an obspy.core.trace.Trace object."
107        raise ValueError(msg)
108
109    if not rtmemory_list:
110        rtmemory_list = [RtMemory()]
111
112    sample = trace.data
113    if np.size(sample) < 1:
114        return sample
115
116    delta_time = trace.stats.delta
117
118    rtmemory = rtmemory_list[0]
119
120    # initialize memory object
121    if not rtmemory.initialized:
122        memory_size_input = 0
123        memory_size_output = 1
124        rtmemory.initialize(sample.dtype, memory_size_input,
125                            memory_size_output, 0, 0)
126
127    sum_ = rtmemory.output[0]
128
129    for i in range(np.size(sample)):
130        sum_ += sample[i] * delta_time
131        sample[i] = sum_
132
133    rtmemory.output[0] = sum_
134
135    return sample
136
137
138def differentiate(trace, rtmemory_list=None):
139    """
140    Apply simple differentiation to array data.
141
142    :type trace: :class:`~obspy.core.trace.Trace`
143    :param trace:  :class:`~obspy.core.trace.Trace` object to append to this
144        RtTrace
145    :type rtmemory_list: list of :class:`~obspy.realtime.rtmemory.RtMemory`,
146        optional
147    :param rtmemory_list: Persistent memory used by this process for specified
148        trace.
149    :rtype: NumPy :class:`numpy.ndarray`
150    :return: Processed trace data from appended Trace object.
151    """
152    if not isinstance(trace, Trace):
153        msg = "trace parameter must be an obspy.core.trace.Trace object."
154        raise ValueError(msg)
155
156    if not rtmemory_list:
157        rtmemory_list = [RtMemory()]
158
159    sample = trace.data
160    if np.size(sample) < 1:
161        return(sample)
162
163    delta_time = trace.stats.delta
164
165    rtmemory = rtmemory_list[0]
166
167    # initialize memory object
168    if not rtmemory.initialized:
169        memory_size_input = 1
170        memory_size_output = 0
171        rtmemory.initialize(sample.dtype, memory_size_input,
172                            memory_size_output, 0, 0)
173        # avoid large diff value for first output sample
174        rtmemory.input[0] = sample[0]
175
176    previous_sample = rtmemory.input[0]
177
178    for i in range(np.size(sample)):
179        diff = (sample[i] - previous_sample) / delta_time
180        previous_sample = sample[i]
181        sample[i] = diff
182
183    rtmemory.input[0] = previous_sample
184
185    return sample
186
187
188def boxcar(trace, width, rtmemory_list=None):
189    """
190    Apply boxcar smoothing to data in array sample.
191
192    :type trace: :class:`~obspy.core.trace.Trace`
193    :param trace:  :class:`~obspy.core.trace.Trace` object to append to this
194        RtTrace
195    :type width: int
196    :param width: Width in number of sample points for filter.
197    :type rtmemory_list: list of :class:`~obspy.realtime.rtmemory.RtMemory`,
198        optional
199    :param rtmemory_list: Persistent memory used by this process for specified
200        trace.
201    :rtype: NumPy :class:`numpy.ndarray`
202    :return: Processed trace data from appended Trace object.
203    """
204    if not isinstance(trace, Trace):
205        msg = "trace parameter must be an obspy.core.trace.Trace object."
206        raise ValueError(msg)
207
208    if not width > 0:
209        msg = "width parameter not specified or < 1."
210        raise ValueError(msg)
211
212    if not rtmemory_list:
213        rtmemory_list = [RtMemory()]
214
215    sample = trace.data
216
217    rtmemory = rtmemory_list[0]
218
219    # initialize memory object
220    if not rtmemory.initialized:
221        memory_size_input = width
222        memory_size_output = 0
223        rtmemory.initialize(sample.dtype, memory_size_input,
224                            memory_size_output, 0, 0)
225
226    # initialize array for time-series results
227    new_sample = np.zeros(np.size(sample), sample.dtype)
228
229    i = 0
230    i1 = i - width
231    i2 = i  # causal boxcar of width width
232    sum_ = 0.0
233    icount = 0
234    for i in range(np.size(sample)):
235        value = 0.0
236        if (icount == 0):  # first pass, accumulate sum
237            for n in range(i1, i2 + 1):
238                if (n < 0):
239                    value = rtmemory.input[width + n]
240                else:
241                    value = sample[n]
242                sum_ += value
243                icount = icount + 1
244        else:  # later passes, update sum
245            if ((i1 - 1) < 0):
246                value = rtmemory.input[width + (i1 - 1)]
247            else:
248                value = sample[(i1 - 1)]
249            sum_ -= value
250            if (i2 < 0):
251                value = rtmemory.input[width + i2]
252            else:
253                value = sample[i2]
254            sum_ += value
255        if (icount > 0):
256            new_sample[i] = (float)(sum_ / float(icount))
257        else:
258            new_sample[i] = 0.0
259        i1 = i1 + 1
260        i2 = i2 + 1
261
262    rtmemory.update_input(sample)
263
264    return new_sample
265
266
267def tauc(trace, width, rtmemory_list=None):
268    """
269    Calculate instantaneous period in a fixed window (Tau_c).
270
271    .. seealso::
272
273        Implements equations 1-3 in [Allen2003]_ except use a fixed width
274        window instead of decay function.
275
276    :type trace: :class:`~obspy.core.trace.Trace`
277    :param trace:  :class:`~obspy.core.trace.Trace` object to append to this
278        RtTrace
279    :type width: int
280    :param width: Width in number of sample points for tauc window.
281    :type rtmemory_list: list of :class:`~obspy.realtime.rtmemory.RtMemory`,
282        optional
283    :param rtmemory_list: Persistent memory used by this process for specified
284        trace.
285    :rtype: NumPy :class:`numpy.ndarray`
286    :return: Processed trace data from appended Trace object.
287    """
288    if not isinstance(trace, Trace):
289        msg = "trace parameter must be an obspy.core.trace.Trace object."
290        raise ValueError(msg)
291
292    if not width > 0:
293        msg = "tauc: width parameter not specified or < 1."
294        raise ValueError(msg)
295
296    if not rtmemory_list:
297        rtmemory_list = [RtMemory(), RtMemory()]
298
299    sample = trace.data
300    delta_time = trace.stats.delta
301
302    rtmemory = rtmemory_list[0]
303    rtmemory_dval = rtmemory_list[1]
304
305    sample_last = 0.0
306
307    # initialize memory object
308    if not rtmemory.initialized:
309        memory_size_input = width
310        memory_size_output = 1
311        rtmemory.initialize(sample.dtype, memory_size_input,
312                            memory_size_output, 0, 0)
313        sample_last = sample[0]
314    else:
315        sample_last = rtmemory.input[width - 1]
316
317    # initialize memory object
318    if not rtmemory_dval.initialized:
319        memory_size_input = width
320        memory_size_output = 1
321        rtmemory_dval.initialize(sample.dtype, memory_size_input,
322                                 memory_size_output, 0, 0)
323
324    new_sample = np.zeros(np.size(sample), sample.dtype)
325    deriv = np.zeros(np.size(sample), sample.dtype)
326
327    # sample_last = rtmemory.input[width - 1]
328    sample_d = 0.0
329    deriv_d = 0.0
330    xval = rtmemory.output[0]
331    dval = rtmemory_dval.output[0]
332
333    for i in range(np.size(sample)):
334
335        sample_d = sample[i]
336        deriv_d = (sample_d - sample_last) / delta_time
337        index_begin = i - width
338        if (index_begin >= 0):
339            xval = xval - (sample[index_begin]) * (sample[index_begin]) \
340                + sample_d * sample_d
341            dval = dval - deriv[index_begin] * deriv[index_begin] \
342                + deriv_d * deriv_d
343        else:
344            index = i
345            xval = xval - rtmemory.input[index] * rtmemory.input[index] \
346                + sample_d * sample_d
347            dval = dval \
348                - rtmemory_dval.input[index] * rtmemory_dval.input[index] \
349                + deriv_d * deriv_d
350        deriv[i] = deriv_d
351        sample_last = sample_d
352        # if (xval > _MIN_FLOAT_VAL &  & dval > _MIN_FLOAT_VAL) {
353        if (dval > _MIN_FLOAT_VAL):
354            new_sample[i] = _TWO_PI * math.sqrt(xval / dval)
355        else:
356            new_sample[i] = 0.0
357
358    # update memory
359    rtmemory.output[0] = xval
360    rtmemory.update_input(sample)
361    rtmemory_dval.output[0] = dval
362    rtmemory_dval.update_input(deriv)
363
364    return new_sample
365
366
367# memory object indices for storing specific values
368_AMP_AT_PICK = 0
369_HAVE_USED_MEMORY = 1
370_FLAG_COMPETE_MWP = 2
371_INT_INT_SUM = 3
372_POLARITY = 4
373_MEMORY_SIZE_OUTPUT = 5
374
375
376def mwpintegral(trace, max_time, ref_time, mem_time=1.0, gain=1.0,
377                rtmemory_list=None):
378    """
379    Calculate Mwp integral on a displacement trace.
380
381    .. seealso:: [Tsuboi1999]_ and [Tsuboi1995]_
382
383    :type trace: :class:`~obspy.core.trace.Trace`
384    :param trace:  :class:`~obspy.core.trace.Trace` object to append to this
385        RtTrace
386    :type max_time: float
387    :param max_time: Maximum time in seconds after ref_time to apply Mwp
388        integration.
389    :type ref_time: :class:`~obspy.core.utcdatetime.UTCDateTime`
390    :param ref_time: Reference date and time of the data sample
391        (e.g. P pick time) at which to begin Mwp integration.
392    :type mem_time: float, optional
393    :param mem_time: Length in seconds of data memory (must be much larger
394        than maximum delay between pick declaration and pick time). Defaults
395        to ``1.0``.
396    :type gain: float, optional
397    :param gain: Nominal gain to convert input displacement trace to meters
398        of ground displacement. Defaults to ``1.0``.
399    :type rtmemory_list: list of :class:`~obspy.realtime.rtmemory.RtMemory`,
400        optional
401    :param rtmemory_list: Persistent memory used by this process for specified
402        trace.
403    :rtype: NumPy :class:`numpy.ndarray`
404    :return: Processed trace data from appended Trace object.
405    """
406    if not isinstance(trace, Trace):
407        msg = "trace parameter must be an obspy.core.trace.Trace object."
408        raise ValueError(msg)
409
410    if not isinstance(ref_time, UTCDateTime):
411        msg = "ref_time must be an obspy.core.utcdatetime.UTCDateTime object."
412        raise ValueError(msg)
413
414    if not max_time >= 0:
415        msg = "max_time parameter not specified or < 0."
416        raise ValueError(msg)
417
418    if not rtmemory_list:
419        rtmemory_list = [RtMemory()]
420
421    sample = trace.data
422    delta_time = trace.stats.delta
423
424    rtmemory = rtmemory_list[0]
425
426    # initialize memory object
427    if not rtmemory.initialized:
428        memory_size_input = int(0.5 + mem_time * trace.stats.sampling_rate)
429        memory_size_output = _MEMORY_SIZE_OUTPUT
430        rtmemory.initialize(sample.dtype, memory_size_input,
431                            memory_size_output, 0, 0)
432
433    new_sample = np.zeros(np.size(sample), sample.dtype)
434
435    ioffset_pick = int(round(
436                       (ref_time - trace.stats.starttime) *
437                       trace.stats.sampling_rate))
438    ioffset_mwp_min = ioffset_pick
439
440    # set reference amplitude
441    if ioffset_mwp_min >= 0 and ioffset_mwp_min < trace.data.size:
442        # value in trace data array
443        rtmemory.output[_AMP_AT_PICK] = trace.data[ioffset_mwp_min]
444    elif ioffset_mwp_min >= -(np.size(rtmemory.input)) and ioffset_mwp_min < 0:
445        # value in memory array
446        index = ioffset_mwp_min + np.size(rtmemory.input)
447        rtmemory.output[_AMP_AT_PICK] = rtmemory.input[index]
448    elif ioffset_mwp_min < -(np.size(rtmemory.input)) \
449            and not rtmemory.output[_HAVE_USED_MEMORY]:
450        msg = "mem_time not large enough to buffer required input data."
451        raise ValueError(msg)
452    if ioffset_mwp_min < 0 and rtmemory.output[_HAVE_USED_MEMORY]:
453        ioffset_mwp_min = 0
454    else:
455        rtmemory.output[_HAVE_USED_MEMORY] = 1
456    # set Mwp end index corresponding to maximum duration
457    mwp_end_index = int(round(max_time / delta_time))
458    ioffset_mwp_max = mwp_end_index + ioffset_pick
459    if ioffset_mwp_max < trace.data.size:
460        rtmemory.output[_FLAG_COMPETE_MWP] = 1  # will complete
461    if ioffset_mwp_max > trace.data.size:
462        ioffset_mwp_max = trace.data.size
463    # apply double integration, check for extrema
464    mwp_amp_at_pick = rtmemory.output[_AMP_AT_PICK]
465    mwp_int_int_sum = rtmemory.output[_INT_INT_SUM]
466    polarity = rtmemory.output[_POLARITY]
467    amplitude = 0.0
468    for n in range(ioffset_mwp_min, ioffset_mwp_max):
469        if n >= 0:
470            amplitude = trace.data[n]
471        elif n >= -(np.size(rtmemory.input)):
472            # value in memory array
473            index = n + np.size(rtmemory.input)
474            amplitude = rtmemory.input[index]
475        else:
476            msg = "Error: Mwp: attempt to access rtmemory.input array of " + \
477                "size=%d at invalid index=%d: this should not happen!" % \
478                (np.size(rtmemory.input), n + np.size(rtmemory.input))
479            print(msg)
480            continue  # should never reach here
481        disp_amp = amplitude - mwp_amp_at_pick
482        # check displacement polarity
483        if disp_amp >= 0.0:  # pos
484            # check if past extremum
485            if polarity < 0:  # passed from neg to pos displacement
486                mwp_int_int_sum *= -1.0
487                mwp_int_int_sum = 0
488            polarity = 1
489        elif disp_amp < 0.0:  # neg
490            # check if past extremum
491            if polarity > 0:  # passed from pos to neg displacement
492                mwp_int_int_sum = 0
493            polarity = -1
494        mwp_int_int_sum += (amplitude - mwp_amp_at_pick) * delta_time / gain
495        new_sample[n] = mwp_int_int_sum
496
497    rtmemory.output[_INT_INT_SUM] = mwp_int_int_sum
498    rtmemory.output[_POLARITY] = polarity
499
500    # update memory
501    rtmemory.update_input(sample)
502
503    return new_sample
504
505
506MWP_INVALID = -9.9
507# 4.213e19 - Tsuboi 1995, 1999
508MWP_CONST = 4.0 * _PI  # 4 PI
509MWP_CONST *= 3400.0  # rho
510MWP_CONST *= 7900.0 * 7900.0 * 7900.0  # Pvel**3
511MWP_CONST *= 2.0  # FP average radiation pattern
512MWP_CONST *= (10000.0 / 90.0)  # distance deg -> km
513MWP_CONST *= 1000.0  # distance km -> meters
514# https://mail.python.org/pipermail/python-list/2010-February/567089.html, ff.
515try:
516    FLOAT_MIN = sys.float_info.min
517except AttributeError:
518    FLOAT_MIN = 1.1e-37
519
520
521def calculate_mwp_mag(peak, epicentral_distance):
522    """
523    Calculate Mwp magnitude.
524
525    .. seealso:: [Tsuboi1999]_ and [Tsuboi1995]_
526
527    :type peak: float
528    :param peak: Peak value of integral of displacement seismogram.
529    :type epicentral_distance: float
530    :param epicentral_distance: Great-circle epicentral distance from station
531        in degrees.
532    :rtype: float
533    :returns: Calculated Mwp magnitude.
534    """
535    moment = MWP_CONST * peak * epicentral_distance
536    mwp_mag = MWP_INVALID
537    if moment > FLOAT_MIN:
538        mwp_mag = (2.0 / 3.0) * (math.log10(moment) - 9.1)
539    return mwp_mag
540
541
542def kurtosis(trace, win=3.0, rtmemory_list=None):
543    """
544    Apply recursive kurtosis calculation on data.
545
546    Recursive kurtosis is computed using the [ChassandeMottin2002]_
547    formulation adjusted to give the kurtosis of a Gaussian distribution = 0.0.
548
549    :type trace: :class:`~obspy.core.trace.Trace`
550    :param trace: :class:`~obspy.core.trace.Trace` object to append to this
551        RtTrace
552    :type win: float, optional
553    :param win: window length in seconds for the kurtosis (default is 3.0 s)
554    :type rtmemory_list: list of :class:`~obspy.realtime.rtmemory.RtMemory`,
555        optional
556    :param rtmemory_list: Persistent memory used by this process for specified
557        trace
558    :rtype: NumPy :class:`numpy.ndarray`
559    :return: Processed trace data from appended Trace object
560    """
561    if not isinstance(trace, Trace):
562        msg = "Trace parameter must be an obspy.core.trace.Trace object."
563        raise ValueError(msg)
564
565    # if this is the first appended trace, the rtmemory_list will be None
566    if not rtmemory_list:
567        rtmemory_list = [RtMemory(), RtMemory(), RtMemory()]
568
569    # deal with case of empty trace
570    sample = trace.data
571    if np.size(sample) < 1:
572        return sample
573
574    # get simple info from trace
575    npts = len(sample)
576    dt = trace.stats.delta
577
578    # set some constants for the kurtosis calculation
579    c_1 = dt / float(win)
580    a1 = 1.0 - c_1
581    c_2 = (1.0 - a1 * a1) / 2.0
582    bias = -3 * c_1 - 3.0
583
584    # prepare the output array
585    kappa4 = np.empty(npts, sample.dtype)
586
587    # initialize the real-time memory needed to store
588    # the recursive kurtosis coefficients until the
589    # next bloc of data is added
590    rtmemory_mu1 = rtmemory_list[0]
591    rtmemory_mu2 = rtmemory_list[1]
592    rtmemory_k4_bar = rtmemory_list[2]
593
594    # there are three memory objects, one for each "last" coefficient
595    # that needs carrying over
596    # initialize mu1_last to 0
597    if not rtmemory_mu1.initialized:
598        memory_size_input = 1
599        memory_size_output = 0
600        rtmemory_mu1.initialize(sample.dtype, memory_size_input,
601                                memory_size_output, 0, 0)
602
603    # initialize mu2_last (sigma) to 1
604    if not rtmemory_mu2.initialized:
605        memory_size_input = 1
606        memory_size_output = 0
607        rtmemory_mu2.initialize(sample.dtype, memory_size_input,
608                                memory_size_output, 1, 0)
609
610    # initialize k4_bar_last to 0
611    if not rtmemory_k4_bar.initialized:
612        memory_size_input = 1
613        memory_size_output = 0
614        rtmemory_k4_bar.initialize(sample.dtype, memory_size_input,
615                                   memory_size_output, 0, 0)
616
617    mu1_last = rtmemory_mu1.input[0]
618    mu2_last = rtmemory_mu2.input[0]
619    k4_bar_last = rtmemory_k4_bar.input[0]
620
621    # do recursive kurtosis
622    for i in range(npts):
623        mu1 = a1 * mu1_last + c_1 * sample[i]
624        dx2 = (sample[i] - mu1_last) * (sample[i] - mu1_last)
625        mu2 = a1 * mu2_last + c_2 * dx2
626        dx2 = dx2 / mu2_last
627        k4_bar = (1 + c_1 - 2 * c_1 * dx2) * k4_bar_last + c_1 * dx2 * dx2
628        kappa4[i] = k4_bar + bias
629        mu1_last = mu1
630        mu2_last = mu2
631        k4_bar_last = k4_bar
632
633    rtmemory_mu1.input[0] = mu1_last
634    rtmemory_mu2.input[0] = mu2_last
635    rtmemory_k4_bar.input[0] = k4_bar_last
636
637    return kappa4
638