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