1#!/usr/bin/env python 2# -*- coding: utf-8 -*- 3""" 4Classes related to instrument responses. 5 6:copyright: 7 Lion Krischer (krischer@geophysik.uni-muenchen.de), 2013 8:license: 9 GNU Lesser General Public License, Version 3 10 (https://www.gnu.org/copyleft/lesser.html) 11""" 12from __future__ import (absolute_import, division, print_function, 13 unicode_literals) 14from future.builtins import * # NOQA 15 16import copy 17import ctypes as C # NOQA 18from collections import defaultdict 19from copy import deepcopy 20import itertools 21from math import pi 22import warnings 23 24import numpy as np 25import scipy.interpolate 26 27from .. import compatibility 28from obspy.core.util.base import ComparingObject 29from obspy.core.util.deprecation_helpers import ObsPyDeprecationWarning 30from obspy.core.util.obspy_types import (ComplexWithUncertainties, 31 FloatWithUncertainties, 32 FloatWithUncertaintiesAndUnit, 33 ObsPyException, 34 ZeroSamplingRate) 35 36from .util import Angle, Frequency 37 38 39class ResponseStage(ComparingObject): 40 """ 41 From the StationXML Definition: 42 This complex type represents channel response and covers SEED 43 blockettes 53 to 56. 44 """ 45 def __init__(self, stage_sequence_number, stage_gain, 46 stage_gain_frequency, input_units, output_units, 47 resource_id=None, resource_id2=None, name=None, 48 input_units_description=None, 49 output_units_description=None, description=None, 50 decimation_input_sample_rate=None, decimation_factor=None, 51 decimation_offset=None, decimation_delay=None, 52 decimation_correction=None): 53 """ 54 :type stage_sequence_number: int 55 :param stage_sequence_number: Stage sequence number, greater or equal 56 to zero. This is used in all the response SEED blockettes. 57 :type stage_gain: float 58 :param stage_gain: Value of stage gain. 59 :type stage_gain_frequency: float 60 :param stage_gain_frequency: Frequency of stage gain. 61 :type input_units: str 62 :param input_units: The units of the data as input from the 63 perspective of data acquisition. After correcting data for this 64 response, these would be the resulting units. 65 Name of units, e.g. "M/S", "V", "PA". 66 :type output_units: str 67 :param output_units: The units of the data as output from the 68 perspective of data acquisition. These would be the units of the 69 data prior to correcting for this response. 70 Name of units, e.g. "M/S", "V", "PA". 71 :type resource_id: str 72 :param resource_id: This field contains a string that should serve as a 73 unique resource identifier. This identifier can be interpreted 74 differently depending on the data center/software that generated 75 the document. Also, we recommend to use something like 76 GENERATOR:Meaningful ID. As a common behavior equipment with the 77 same ID should contains the same information/be derived from the 78 same base instruments. 79 :type resource_id2: str 80 :param resource_id2: This field contains a string that should serve as 81 a unique resource identifier. Resource identifier of the subgroup 82 of the response stage that varies across different response stage 83 types (e.g. the poles and zeros part or the FIR part). 84 :type name: str 85 :param name: A name given to the filter stage. 86 :type input_units_description: str, optional 87 :param input_units_description: The units of the data as input from the 88 perspective of data acquisition. After correcting data for this 89 response, these would be the resulting units. 90 Description of units, e.g. "Velocity in meters per second", 91 "Volts", "Pascals". 92 :type output_units_description: str, optional 93 :param output_units_description: The units of the data as output from 94 the perspective of data acquisition. These would be the units of 95 the data prior to correcting for this response. 96 Description of units, e.g. "Velocity in meters per second", 97 "Volts", "Pascals". 98 :type description: str, optional 99 :param description: A short description of of the filter. 100 :type decimation_input_sample_rate: float, optional 101 :param decimation_input_sample_rate: The sampling rate before the 102 decimation in samples per second. 103 :type decimation_factor: int, optional 104 :param decimation_factor: The applied decimation factor. 105 :type decimation_offset: int, optional 106 :param decimation_offset: The sample chosen for use. 0 denotes the 107 first sample, 1 the second, and so forth. 108 :type decimation_delay: float, optional 109 :param decimation_delay: The estimated pure delay from the decimation. 110 :type decimation_correction: float, optional 111 :param decimation_correction: The time shift applied to correct for the 112 delay at this stage. 113 114 .. note:: 115 The stage gain (or stage sensitivity) is the gain at the stage of 116 the encapsulating response element and corresponds to SEED 117 blockette 58. In the SEED convention, stage 0 gain represents the 118 overall sensitivity of the channel. In this schema, stage 0 gains 119 are allowed but are considered deprecated. Overall sensitivity 120 should be specified in the InstrumentSensitivity element. 121 """ 122 self.stage_sequence_number = stage_sequence_number 123 self.input_units = input_units 124 self.output_units = output_units 125 self.input_units_description = input_units_description 126 self.output_units_description = output_units_description 127 self.resource_id = resource_id 128 self.resource_id2 = resource_id2 129 self.stage_gain = stage_gain 130 self.stage_gain_frequency = stage_gain_frequency 131 self.name = name 132 self.description = description 133 self.decimation_input_sample_rate = \ 134 Frequency(decimation_input_sample_rate) \ 135 if decimation_input_sample_rate is not None else None 136 self.decimation_factor = decimation_factor 137 self.decimation_offset = decimation_offset 138 self.decimation_delay = \ 139 FloatWithUncertaintiesAndUnit(decimation_delay) \ 140 if decimation_delay is not None else None 141 self.decimation_correction = \ 142 FloatWithUncertaintiesAndUnit(decimation_correction) \ 143 if decimation_correction is not None else None 144 145 def __str__(self): 146 ret = ( 147 "Response type: {response_type}, Stage Sequence Number: " 148 "{response_stage}\n" 149 "{name_desc}" 150 "{resource_id}" 151 "\tFrom {input_units}{input_desc} to {output_units}{output_desc}\n" 152 "\tStage gain: {gain}, defined at {gain_freq} Hz\n" 153 "{decimation}").format( 154 response_type=self.__class__.__name__, 155 response_stage=self.stage_sequence_number, 156 gain=self.stage_gain if self.stage_gain is not None else "UNKNOWN", 157 gain_freq=("%.2f" % self.stage_gain_frequency) if 158 self.stage_gain_frequency is not None else "UNKNOWN", 159 name_desc="\t%s %s\n" % ( 160 self.name, "(%s)" % self.description 161 if self.description else "") if self.name else "", 162 resource_id=("\tResource Id: %s" % self.resource_id 163 if self.resource_id else ""), 164 input_units=self.input_units if self.input_units else "UNKNOWN", 165 input_desc=(" (%s)" % self.input_units_description 166 if self.input_units_description else ""), 167 output_units=self.output_units if self.output_units else "UNKNOWN", 168 output_desc=(" (%s)" % self.output_units_description 169 if self.output_units_description else ""), 170 decimation=( 171 "\tDecimation:\n\t\tInput Sample Rate: %.2f Hz\n\t\t" 172 "Decimation Factor: %i\n\t\tDecimation Offset: %i\n\t\t" 173 "Decimation Delay: %.2f\n\t\tDecimation Correction: %.2f" % ( 174 self.decimation_input_sample_rate, self.decimation_factor, 175 self.decimation_offset, self.decimation_delay, 176 self.decimation_correction) 177 if self.decimation_input_sample_rate is not None else "")) 178 return ret.strip() 179 180 def _repr_pretty_(self, p, cycle): 181 p.text(str(self)) 182 183 184class PolesZerosResponseStage(ResponseStage): 185 """ 186 From the StationXML Definition: 187 Response: complex poles and zeros. Corresponds to SEED blockette 53. 188 189 The response stage is used for the analog stages of the filter system and 190 the description of infinite impulse response (IIR) digital filters. 191 192 Has all the arguments of the parent class 193 :class:`~obspy.core.inventory.response.ResponseStage` and the following: 194 195 :type pz_transfer_function_type: str 196 :param pz_transfer_function_type: A string describing the type of transfer 197 function. Can be one of: 198 199 * ``LAPLACE (RADIANS/SECOND)`` 200 * ``LAPLACE (HERTZ)`` 201 * ``DIGITAL (Z-TRANSFORM)`` 202 203 The function tries to match inputs to one of three types if it can. 204 :type normalization_frequency: float 205 :param normalization_frequency: The frequency at which the normalization 206 factor is normalized. 207 :type zeros: list of complex 208 :param zeros: All zeros of the stage. 209 :type poles: list of complex 210 :param poles: All poles of the stage. 211 :type normalization_factor: float, optional 212 :param normalization_factor: 213 """ 214 def __init__(self, stage_sequence_number, stage_gain, 215 stage_gain_frequency, input_units, output_units, 216 pz_transfer_function_type, 217 normalization_frequency, zeros, poles, 218 normalization_factor=1.0, resource_id=None, resource_id2=None, 219 name=None, input_units_description=None, 220 output_units_description=None, description=None, 221 decimation_input_sample_rate=None, decimation_factor=None, 222 decimation_offset=None, decimation_delay=None, 223 decimation_correction=None): 224 # Set the Poles and Zeros specific attributes. Special cases are 225 # handled by properties. 226 self.pz_transfer_function_type = pz_transfer_function_type 227 self.normalization_frequency = normalization_frequency 228 self.normalization_factor = float(normalization_factor) 229 self.zeros = zeros 230 self.poles = poles 231 super(PolesZerosResponseStage, self).__init__( 232 stage_sequence_number=stage_sequence_number, 233 input_units=input_units, 234 output_units=output_units, 235 input_units_description=input_units_description, 236 output_units_description=output_units_description, 237 resource_id=resource_id, resource_id2=resource_id2, 238 stage_gain=stage_gain, 239 stage_gain_frequency=stage_gain_frequency, name=name, 240 description=description, 241 decimation_input_sample_rate=decimation_input_sample_rate, 242 decimation_factor=decimation_factor, 243 decimation_offset=decimation_offset, 244 decimation_delay=decimation_delay, 245 decimation_correction=decimation_correction) 246 247 def __str__(self): 248 ret = super(PolesZerosResponseStage, self).__str__() 249 ret += ( 250 "\n" 251 "\tTransfer function type: {transfer_fct_type}\n" 252 "\tNormalization factor: {norm_fact:g}, " 253 "Normalization frequency: {norm_freq:.2f} Hz\n" 254 "\tPoles: {poles}\n" 255 "\tZeros: {zeros}").format( 256 transfer_fct_type=self.pz_transfer_function_type, 257 norm_fact=self.normalization_factor, 258 norm_freq=self.normalization_frequency, 259 poles=", ".join(map(str, self.poles)), 260 zeros=", ".join(map(str, self.zeros))) 261 return ret 262 263 def _repr_pretty_(self, p, cycle): 264 p.text(str(self)) 265 266 @property 267 def zeros(self): 268 return self._zeros 269 270 @zeros.setter 271 def zeros(self, value): 272 value = list(value) 273 for i, x in enumerate(value): 274 if not isinstance(x, ComplexWithUncertainties): 275 value[i] = ComplexWithUncertainties(x) 276 self._zeros = value 277 278 @property 279 def poles(self): 280 return self._poles 281 282 @poles.setter 283 def poles(self, value): 284 value = list(value) 285 for i, x in enumerate(value): 286 if not isinstance(x, ComplexWithUncertainties): 287 value[i] = ComplexWithUncertainties(x) 288 self._poles = value 289 290 @property 291 def normalization_frequency(self): 292 return self._normalization_frequency 293 294 @normalization_frequency.setter 295 def normalization_frequency(self, value): 296 self._normalization_frequency = Frequency(value) 297 298 @property 299 def pz_transfer_function_type(self): 300 return self._pz_transfer_function_type 301 302 @pz_transfer_function_type.setter 303 def pz_transfer_function_type(self, value): 304 """ 305 Setter for the transfer function type. 306 307 Rather permissive but should make it less awkward to use. 308 """ 309 msg = ("'%s' is not a valid value for 'pz_transfer_function_type'. " 310 "Valid one are:\n" 311 "\tLAPLACE (RADIANS/SECOND)\n" 312 "\tLAPLACE (HERTZ)\n" 313 "\tDIGITAL (Z-TRANSFORM)") % value 314 value = value.lower() 315 if "laplace" in value: 316 if "radian" in value: 317 self._pz_transfer_function_type = "LAPLACE (RADIANS/SECOND)" 318 elif "hertz" in value or "hz" in value: 319 self._pz_transfer_function_type = "LAPLACE (HERTZ)" 320 else: 321 raise ValueError(msg) 322 elif "digital" in value: 323 self._pz_transfer_function_type = "DIGITAL (Z-TRANSFORM)" 324 else: 325 raise ValueError(msg) 326 327 328class CoefficientsTypeResponseStage(ResponseStage): 329 """ 330 This response type can describe coefficients for FIR filters. Laplace 331 transforms and IIR filters can also be expressed using this type but should 332 rather be described using the PolesZerosResponseStage class. Effectively 333 corresponds to SEED blockette 54. 334 335 Has all the arguments of the parent class 336 :class:`~obspy.core.inventory.response.ResponseStage` and the following: 337 338 :type cf_transfer_function_type: str 339 :param cf_transfer_function_type: A string describing the type of transfer 340 function. Can be one of: 341 342 * ``ANALOG (RADIANS/SECOND)`` 343 * ``ANALOG (HERTZ)`` 344 * ``DIGITAL`` 345 346 The function tries to match inputs to one of three types if it can. 347 :type numerator: list of 348 :class:`~obspy.core.util.obspy_types.CoefficientWithUncertainties` 349 :param numerator: Numerator of the coefficient response stage. 350 :type denominator: list of 351 :class:`~obspy.core.util.obspy_types.CoefficientWithUncertainties` 352 :param denominator: Denominator of the coefficient response stage. 353 """ 354 def __init__(self, stage_sequence_number, stage_gain, 355 stage_gain_frequency, input_units, output_units, 356 cf_transfer_function_type, resource_id=None, 357 resource_id2=None, name=None, numerator=None, 358 denominator=None, input_units_description=None, 359 output_units_description=None, description=None, 360 decimation_input_sample_rate=None, decimation_factor=None, 361 decimation_offset=None, decimation_delay=None, 362 decimation_correction=None): 363 # Set the Coefficients type specific attributes. Special cases are 364 # handled by properties. 365 self.cf_transfer_function_type = cf_transfer_function_type 366 self.numerator = numerator 367 self.denominator = denominator 368 super(CoefficientsTypeResponseStage, self).__init__( 369 stage_sequence_number=stage_sequence_number, 370 input_units=input_units, 371 output_units=output_units, 372 input_units_description=input_units_description, 373 output_units_description=output_units_description, 374 resource_id=resource_id, resource_id2=resource_id2, 375 stage_gain=stage_gain, 376 stage_gain_frequency=stage_gain_frequency, name=name, 377 description=description, 378 decimation_input_sample_rate=decimation_input_sample_rate, 379 decimation_factor=decimation_factor, 380 decimation_offset=decimation_offset, 381 decimation_delay=decimation_delay, 382 decimation_correction=decimation_correction) 383 384 def __str__(self): 385 ret = super(CoefficientsTypeResponseStage, self).__str__() 386 ret += ( 387 "\n" 388 "\tTransfer function type: {transfer_fct_type}\n" 389 "\tContains {num_count} numerators and {den_count} denominators")\ 390 .format( 391 transfer_fct_type=self.cf_transfer_function_type, 392 num_count=len(self.numerator), den_count=len(self.denominator)) 393 return ret 394 395 def _repr_pretty_(self, p, cycle): 396 p.text(str(self)) 397 398 @property 399 def numerator(self): 400 return self._numerator 401 402 @numerator.setter 403 def numerator(self, value): 404 if value == []: 405 self._numerator = [] 406 return 407 value = list(value) if isinstance( 408 value, compatibility.collections_abc.Iterable) else [value] 409 if any(getattr(x, 'unit', None) is not None for x in value): 410 msg = 'Setting Numerator/Denominator with a unit is deprecated.' 411 warnings.warn(msg, ObsPyDeprecationWarning) 412 for _i, x in enumerate(value): 413 if not isinstance(x, CoefficientWithUncertainties): 414 value[_i] = CoefficientWithUncertainties(x) 415 self._numerator = value 416 417 @property 418 def denominator(self): 419 return self._denominator 420 421 @denominator.setter 422 def denominator(self, value): 423 if value == []: 424 self._denominator = [] 425 return 426 value = list(value) if isinstance( 427 value, compatibility.collections_abc.Iterable) else [value] 428 if any(getattr(x, 'unit', None) is not None for x in value): 429 msg = 'Setting Numerator/Denominator with a unit is deprecated.' 430 warnings.warn(msg, ObsPyDeprecationWarning) 431 for _i, x in enumerate(value): 432 if not isinstance(x, CoefficientWithUncertainties): 433 value[_i] = CoefficientWithUncertainties(x) 434 self._denominator = value 435 436 @property 437 def cf_transfer_function_type(self): 438 return self._cf_transfer_function_type 439 440 @cf_transfer_function_type.setter 441 def cf_transfer_function_type(self, value): 442 """ 443 Setter for the transfer function type. 444 445 Rather permissive but should make it less awkward to use. 446 """ 447 msg = ("'%s' is not a valid value for 'cf_transfer_function_type'. " 448 "Valid one are:\n" 449 "\tANALOG (RADIANS/SECOND)\n" 450 "\tANALOG (HERTZ)\n" 451 "\tDIGITAL") % value 452 value = value.lower() 453 if "analog" in value: 454 if "rad" in value: 455 self._cf_transfer_function_type = "ANALOG (RADIANS/SECOND)" 456 elif "hertz" in value or "hz" in value: 457 self._cf_transfer_function_type = "ANALOG (HERTZ)" 458 else: 459 raise ValueError(msg) 460 elif "digital" in value: 461 self._cf_transfer_function_type = "DIGITAL" 462 else: 463 raise ValueError(msg) 464 465 466class ResponseListResponseStage(ResponseStage): 467 """ 468 This response type gives a list of frequency, amplitude and phase value 469 pairs. Effectively corresponds to SEED blockette 55. 470 471 Has all the arguments of the parent class 472 :class:`~obspy.core.inventory.response.ResponseStage` and the following: 473 474 :type response_list_elements: list of 475 :class:`~obspy.core.inventory.response.ResponseListElement` 476 :param response_list_elements: A list of single discrete frequency, 477 amplitude and phase response values. 478 """ 479 def __init__(self, stage_sequence_number, stage_gain, 480 stage_gain_frequency, input_units, output_units, 481 resource_id=None, resource_id2=None, name=None, 482 response_list_elements=None, 483 input_units_description=None, output_units_description=None, 484 description=None, decimation_input_sample_rate=None, 485 decimation_factor=None, decimation_offset=None, 486 decimation_delay=None, decimation_correction=None): 487 self.response_list_elements = response_list_elements or [] 488 super(ResponseListResponseStage, self).__init__( 489 stage_sequence_number=stage_sequence_number, 490 input_units=input_units, 491 output_units=output_units, 492 input_units_description=input_units_description, 493 output_units_description=output_units_description, 494 resource_id=resource_id, resource_id2=resource_id2, 495 stage_gain=stage_gain, 496 stage_gain_frequency=stage_gain_frequency, name=name, 497 description=description, 498 decimation_input_sample_rate=decimation_input_sample_rate, 499 decimation_factor=decimation_factor, 500 decimation_offset=decimation_offset, 501 decimation_delay=decimation_delay, 502 decimation_correction=decimation_correction) 503 504 505class ResponseListElement(ComparingObject): 506 """ 507 Describes the amplitude and phase response value for a discrete frequency 508 value. 509 """ 510 def __init__(self, frequency, amplitude, phase): 511 """ 512 :type frequency: float 513 :param frequency: The frequency for which the response is valid. 514 :type amplitude: float 515 :param amplitude: The value for the amplitude response at this 516 frequency. 517 :type phase: float 518 :param phase: The value for the phase response at this frequency. 519 """ 520 self.frequency = frequency 521 self.amplitude = amplitude 522 self.phase = phase 523 524 @property 525 def frequency(self): 526 return self._frequency 527 528 @frequency.setter 529 def frequency(self, value): 530 if not isinstance(value, Frequency): 531 value = Frequency(value) 532 self._frequency = value 533 534 @property 535 def amplitude(self): 536 return self._amplitude 537 538 @amplitude.setter 539 def amplitude(self, value): 540 if not isinstance(value, FloatWithUncertaintiesAndUnit): 541 value = FloatWithUncertaintiesAndUnit(value) 542 self._amplitude = value 543 544 @property 545 def phase(self): 546 return self._phase 547 548 @phase.setter 549 def phase(self, value): 550 if not isinstance(value, Angle): 551 value = Angle(value) 552 self._phase = value 553 554 555class FIRResponseStage(ResponseStage): 556 """ 557 From the StationXML Definition: 558 Response: FIR filter. Corresponds to SEED blockette 61. FIR filters are 559 also commonly documented using the CoefficientsType element. 560 561 Has all the arguments of the parent class 562 :class:`~obspy.core.inventory.response.ResponseStage` and the following: 563 564 :type symmetry: str 565 :param symmetry: A string describing the symmetry. Can be one of: 566 567 * ``NONE`` 568 * ``EVEN`` 569 * ``ODD`` 570 571 :type coefficients: list of floats 572 :param coefficients: List of FIR coefficients. 573 """ 574 def __init__(self, stage_sequence_number, stage_gain, 575 stage_gain_frequency, input_units, output_units, 576 symmetry="NONE", resource_id=None, resource_id2=None, 577 name=None, 578 coefficients=None, input_units_description=None, 579 output_units_description=None, description=None, 580 decimation_input_sample_rate=None, decimation_factor=None, 581 decimation_offset=None, decimation_delay=None, 582 decimation_correction=None): 583 self._symmetry = symmetry 584 self.coefficients = coefficients or [] 585 super(FIRResponseStage, self).__init__( 586 stage_sequence_number=stage_sequence_number, 587 input_units=input_units, 588 output_units=output_units, 589 input_units_description=input_units_description, 590 output_units_description=output_units_description, 591 resource_id=resource_id, resource_id2=resource_id2, 592 stage_gain=stage_gain, 593 stage_gain_frequency=stage_gain_frequency, name=name, 594 description=description, 595 decimation_input_sample_rate=decimation_input_sample_rate, 596 decimation_factor=decimation_factor, 597 decimation_offset=decimation_offset, 598 decimation_delay=decimation_delay, 599 decimation_correction=decimation_correction) 600 601 @property 602 def symmetry(self): 603 return self._symmetry 604 605 @symmetry.setter 606 def symmetry(self, value): 607 value = str(value).upper() 608 allowed = ("NONE", "EVEN", "ODD") 609 if value not in allowed: 610 msg = ("Value '%s' for FIR Response symmetry not allowed. " 611 "Possible values are: '%s'") 612 msg = msg % (value, "', '".join(allowed)) 613 raise ValueError(msg) 614 self._symmetry = value 615 616 @property 617 def coefficients(self): 618 return self._coefficients 619 620 @coefficients.setter 621 def coefficients(self, value): 622 new_values = [] 623 for x in value: 624 if not isinstance(x, FilterCoefficient): 625 x = FilterCoefficient(x) 626 new_values.append(x) 627 self._coefficients = new_values 628 629 630class PolynomialResponseStage(ResponseStage): 631 """ 632 From the StationXML Definition: 633 Response: expressed as a polynomial (allows non-linear sensors to be 634 described). Corresponds to SEED blockette 62. Can be used to describe a 635 stage of acquisition or a complete system. 636 637 Has all the arguments of the parent class 638 :class:`~obspy.core.inventory.response.ResponseStage` and the following: 639 640 :type approximation_type: str 641 :param approximation_type: Approximation type. Currently restricted to 642 'MACLAURIN' by StationXML definition. 643 :type frequency_lower_bound: float 644 :param frequency_lower_bound: Lower frequency bound. 645 :type frequency_upper_bound: float 646 :param frequency_upper_bound: Upper frequency bound. 647 :type approximation_lower_bound: float 648 :param approximation_lower_bound: Lower bound of approximation. 649 :type approximation_upper_bound: float 650 :param approximation_upper_bound: Upper bound of approximation. 651 :type maximum_error: float 652 :param maximum_error: Maximum error. 653 :type coefficients: list of floats 654 :param coefficients: List of polynomial coefficients. 655 """ 656 def __init__(self, stage_sequence_number, stage_gain, 657 stage_gain_frequency, input_units, output_units, 658 frequency_lower_bound, 659 frequency_upper_bound, approximation_lower_bound, 660 approximation_upper_bound, maximum_error, coefficients, 661 approximation_type='MACLAURIN', resource_id=None, 662 resource_id2=None, name=None, 663 input_units_description=None, 664 output_units_description=None, description=None, 665 decimation_input_sample_rate=None, decimation_factor=None, 666 decimation_offset=None, decimation_delay=None, 667 decimation_correction=None): 668 # XXX remove stage_gain and stage_gain_frequency completely, since 669 # changes in StationXML 1.1? 670 self._approximation_type = approximation_type 671 self.frequency_lower_bound = frequency_lower_bound 672 self.frequency_upper_bound = frequency_upper_bound 673 self.approximation_lower_bound = approximation_lower_bound 674 self.approximation_upper_bound = approximation_upper_bound 675 self.maximum_error = maximum_error 676 self.coefficients = coefficients 677 # XXX StationXML 1.1 does not allow stage gain in Polynomial response 678 # stages. Maybe we should we warn here.. but this could get very 679 # verbose when reading StationXML 1.0 files, so maybe not 680 super(PolynomialResponseStage, self).__init__( 681 stage_sequence_number=stage_sequence_number, 682 input_units=input_units, 683 output_units=output_units, 684 input_units_description=input_units_description, 685 output_units_description=output_units_description, 686 resource_id=resource_id, resource_id2=resource_id2, 687 stage_gain=stage_gain, 688 stage_gain_frequency=stage_gain_frequency, name=name, 689 description=description, 690 decimation_input_sample_rate=decimation_input_sample_rate, 691 decimation_factor=decimation_factor, 692 decimation_offset=decimation_offset, 693 decimation_delay=decimation_delay, 694 decimation_correction=decimation_correction) 695 696 @property 697 def approximation_type(self): 698 return self._approximation_type 699 700 @approximation_type.setter 701 def approximation_type(self, value): 702 value = str(value).upper() 703 allowed = ("MACLAURIN",) 704 if value not in allowed: 705 msg = ("Value '%s' for polynomial response approximation type not " 706 "allowed. Possible values are: '%s'") 707 msg = msg % (value, "', '".join(allowed)) 708 raise ValueError(msg) 709 self._approximation_type = value 710 711 @property 712 def coefficients(self): 713 return self._coefficients 714 715 @coefficients.setter 716 def coefficients(self, value): 717 new_values = [] 718 for x in value: 719 if not isinstance(x, CoefficientWithUncertainties): 720 x = CoefficientWithUncertainties(x) 721 new_values.append(x) 722 self._coefficients = new_values 723 724 def __str__(self): 725 ret = super(PolynomialResponseStage, self).__str__() 726 ret += ( 727 "\n" 728 "\tPolynomial approximation type: {approximation_type}\n" 729 "\tFrequency lower bound: {lower_freq_bound}\n" 730 "\tFrequency upper bound: {upper_freq_bound}\n" 731 "\tApproximation lower bound: {lower_approx_bound}\n" 732 "\tApproximation upper bound: {upper_approx_bound}\n" 733 "\tMaximum error: {max_error}\n" 734 "\tNumber of coefficients: {coeff_count}".format( 735 approximation_type=self._approximation_type, 736 lower_freq_bound=self.frequency_lower_bound, 737 upper_freq_bound=self.frequency_upper_bound, 738 lower_approx_bound=self.approximation_lower_bound, 739 upper_approx_bound=self.approximation_upper_bound, 740 max_error=self.maximum_error, 741 coeff_count=len(self.coefficients))) 742 return ret 743 744 745class Response(ComparingObject): 746 """ 747 The root response object. 748 """ 749 def __init__(self, resource_id=None, instrument_sensitivity=None, 750 instrument_polynomial=None, response_stages=None): 751 """ 752 :type resource_id: str 753 :param resource_id: This field contains a string that should serve as a 754 unique resource identifier. This identifier can be interpreted 755 differently depending on the data center/software that generated 756 the document. Also, we recommend to use something like 757 GENERATOR:Meaningful ID. As a common behavior equipment with the 758 same ID should contains the same information/be derived from the 759 same base instruments. 760 :type instrument_sensitivity: 761 :class:`~obspy.core.inventory.response.InstrumentSensitivity` 762 :param instrument_sensitivity: The total sensitivity for the given 763 channel, representing the complete acquisition system expressed as 764 a scalar. 765 :type instrument_polynomial: 766 :class:`~obspy.core.inventory.response.InstrumentPolynomial` 767 :param instrument_polynomial: The total sensitivity for the given 768 channel, representing the complete acquisition system expressed as 769 a polynomial. 770 :type response_stages: list of 771 :class:`~obspy.core.inventory.response.ResponseStage` objects 772 :param response_stages: A list of the response stages. Covers SEED 773 blockettes 53 to 56. 774 """ 775 self.resource_id = resource_id 776 self.instrument_sensitivity = instrument_sensitivity 777 self.instrument_polynomial = instrument_polynomial 778 if response_stages is None: 779 self.response_stages = [] 780 elif hasattr(response_stages, "__iter__"): 781 self.response_stages = response_stages 782 else: 783 msg = "response_stages must be an iterable." 784 raise ValueError(msg) 785 786 def _attempt_to_fix_units(self): 787 """ 788 Internal helper function that will add units to gain only stages based 789 on the units of surrounding stages. 790 791 Should be called when parsing from file formats that don't have units 792 for identity stages. 793 """ 794 previous_output_units = None 795 previous_output_units_description = None 796 797 # Potentially set the input units of the first stage to the units of 798 # the overall sensitivity and the output units of the second stage. 799 if self.response_stages and self.response_stages[0] and \ 800 hasattr(self, "instrument_sensitivity"): 801 s = self.instrument_sensitivity 802 803 if s: 804 if self.response_stages[0].input_units is None: 805 self.response_stages[0].input_units = s.input_units 806 if self.response_stages[0].input_units_description is None: 807 self.response_stages[0].input_units_description = \ 808 s.input_units_description 809 810 if len(self.response_stages) >= 2 and self.response_stages[1]: 811 if self.response_stages[0].output_units is None: 812 self.response_stages[0].output_units = \ 813 self.response_stages[1].input_units 814 if self.response_stages[0].output_units_description is None: 815 self.response_stages[0].output_units_description = \ 816 self.response_stages[1].input_units_description 817 818 # Front to back. 819 for r in self.response_stages: 820 # Only for identity/stage only. 821 if type(r) is ResponseStage: 822 if not r.input_units and not r.output_units and \ 823 previous_output_units: 824 r.input_units = previous_output_units 825 r.output_units = previous_output_units 826 if not r.input_units_description and \ 827 not r.output_units_description \ 828 and previous_output_units_description: 829 r.input_units_description = \ 830 previous_output_units_description 831 r.output_units_description = \ 832 previous_output_units_description 833 834 previous_output_units = r.output_units 835 previous_output_units_description = r.output_units_description 836 837 # Back to front. 838 previous_input_units = None 839 previous_input_units_description = None 840 for r in reversed(self.response_stages): 841 # Only for identity/stage only. 842 if type(r) is ResponseStage: 843 if not r.input_units and not r.output_units and \ 844 previous_input_units: 845 r.input_units = previous_input_units 846 r.output_units = previous_input_units 847 if not r.input_units_description and \ 848 not r.output_units_description \ 849 and previous_input_units_description: 850 r.input_units_description = \ 851 previous_input_units_description 852 r.output_units_description = \ 853 previous_input_units_description 854 855 previous_input_units = r.input_units 856 previous_input_units_description = r.input_units_description 857 858 def get_sampling_rates(self): 859 """ 860 Computes the input and output sampling rates of each stage. 861 862 For well defined files this will just read the decimation attributes 863 of each stage. For others it will attempt to infer missing values 864 from the surrounding stages. 865 866 :returns: A nested dictionary detailing the sampling rates of each 867 response stage. 868 :rtype: dict 869 870 >>> import obspy 871 >>> inv = obspy.read_inventory("AU.MEEK.xml") # doctest: +SKIP 872 >>> inv[0][0][0].response.get_sampling_rates() # doctest: +SKIP 873 {1: {'decimation_factor': 1, 874 'input_sampling_rate': 600.0, 875 'output_sampling_rate': 600.0}, 876 2: {'decimation_factor': 1, 877 'input_sampling_rate': 600.0, 878 'output_sampling_rate': 600.0}, 879 3: {'decimation_factor': 1, 880 'input_sampling_rate': 600.0, 881 'output_sampling_rate': 600.0}, 882 4: {'decimation_factor': 3, 883 'input_sampling_rate': 600.0, 884 'output_sampling_rate': 200.0}, 885 5: {'decimation_factor': 10, 886 'input_sampling_rate': 200.0, 887 'output_sampling_rate': 20.0}} 888 """ 889 # Get all stages, but skip stage 0. 890 stages = [_i.stage_sequence_number for _i in self.response_stages 891 if _i.stage_sequence_number] 892 if not stages: 893 return {} 894 895 if list(range(1, len(stages) + 1)) != stages: 896 raise ValueError("Can only determine sampling rates if response " 897 "stages are in order.") 898 899 # First fill in all the set values. 900 sampling_rates = {} 901 for stage in self.response_stages: 902 input_sr = None 903 output_sr = None 904 factor = None 905 if stage.decimation_input_sample_rate: 906 input_sr = stage.decimation_input_sample_rate 907 if stage.decimation_factor: 908 factor = stage.decimation_factor 909 output_sr = input_sr / float(factor) 910 sampling_rates[stage.stage_sequence_number] = { 911 "input_sampling_rate": input_sr, 912 "output_sampling_rate": output_sr, 913 "decimation_factor": factor} 914 915 # Nothing might be set - just return in that case. 916 if set(itertools.chain.from_iterable(v.values() 917 for v in sampling_rates.values())) == {None}: 918 return sampling_rates 919 920 # Find the first set input sampling rate. The output sampling rate 921 # cannot be set without it. Set all prior input and output sampling 922 # rates to it. 923 for i in stages: 924 sr = sampling_rates[i]["input_sampling_rate"] 925 if sr: 926 for j in range(1, i): 927 sampling_rates[j]["input_sampling_rate"] = sr 928 sampling_rates[j]["output_sampling_rate"] = sr 929 sampling_rates[j]["decimation_factor"] = 1 930 break 931 932 # This should guarantee that the input and output sampling rate of the 933 # the first stage are set. 934 output_sr = sampling_rates[1]["output_sampling_rate"] 935 if not output_sr: # pragma: no cover 936 raise NotImplementedError 937 938 for i in stages: 939 si = sampling_rates[i] 940 if not si["input_sampling_rate"]: 941 si["input_sampling_rate"] = output_sr 942 if not si["output_sampling_rate"]: 943 if not si["decimation_factor"]: 944 si["output_sampling_rate"] = si["input_sampling_rate"] 945 si["decimation_factor"] = 1 946 else: 947 si["output_sampling_rate"] = si["input_sampling_rate"] / \ 948 float(si["decimation_factor"]) 949 if not si["decimation_factor"]: 950 si["decimation_factor"] = int(round( 951 si["input_sampling_rate"] / si["output_sampling_rate"])) 952 output_sr = si["output_sampling_rate"] 953 954 def is_close(a, b): 955 return abs(a - b) < 1e-5 956 957 # Final consistency checks. 958 sr = sampling_rates[stages[0]]["input_sampling_rate"] 959 for i in stages: 960 si = sampling_rates[i] 961 if not is_close(si["input_sampling_rate"], sr): # pragma: no cover 962 msg = ("Input sampling rate of stage %i is inconsistent " 963 "with the previous stages' output sampling rate") 964 warnings.warn(msg % i) 965 966 if not is_close( 967 si["input_sampling_rate"] / si["output_sampling_rate"], 968 si["decimation_factor"]): # pragma: no cover 969 msg = ("Decimation factor in stage %i is inconsistent with " 970 "input and output sampling rates.") 971 warnings.warn(msg % i) 972 sr = si["output_sampling_rate"] 973 974 return sampling_rates 975 976 def recalculate_overall_sensitivity(self, frequency=None): 977 """ 978 Recalculates the overall sensitivity. 979 980 :param frequency: Choose frequency at which to calculate the 981 sensitivity. If not given it will be chosen automatically. 982 """ 983 if not hasattr(self, "instrument_sensitivity"): 984 msg = "Could not find an instrument sensitivity - will not " \ 985 "recalculate the overall sensitivity." 986 raise ValueError(msg) 987 988 if not self.instrument_sensitivity.input_units: 989 msg = "Could not determine input units - will not " \ 990 "recalculate the overall sensitivity." 991 raise ValueError(msg) 992 993 i_u = self.instrument_sensitivity.input_units 994 995 unit_map = { 996 "DISP": ["M"], 997 "VEL": ["M/S", "M/SEC"], 998 "ACC": ["M/S**2", "M/(S**2)", "M/SEC**2", "M/(SEC**2)", 999 "M/S/S"]} 1000 unit = None 1001 for key, value in unit_map.items(): 1002 if i_u and i_u.upper() in value: 1003 unit = key 1004 if not unit: 1005 msg = ("ObsPy does not know how to map unit '%s' to " 1006 "displacement, velocity, or acceleration - overall " 1007 "sensitivity will not be recalculated.") % i_u 1008 raise ValueError(msg) 1009 1010 # Determine frequency if not given. 1011 if frequency is None: 1012 # lookup normalization frequency of sensor's first stage it should 1013 # be in the flat part of the response 1014 stage_one = self.response_stages[0] 1015 try: 1016 frequency = stage_one.normalization_frequency 1017 except AttributeError: 1018 pass 1019 for stage in self.response_stages[::-1]: 1020 # determine sampling rate 1021 try: 1022 sampling_rate = (stage.decimation_input_sample_rate / 1023 stage.decimation_factor) 1024 break 1025 except Exception: 1026 continue 1027 else: 1028 sampling_rate = None 1029 if sampling_rate: 1030 # if sensor's normalization frequency is above 0.5 * nyquist, 1031 # use that instead (e.g. to avoid computing an overall 1032 # sensitivity above nyquist) 1033 nyquist = sampling_rate / 2.0 1034 if frequency: 1035 frequency = min(frequency, nyquist / 2.0) 1036 else: 1037 frequency = nyquist / 2.0 1038 1039 if frequency is None: 1040 msg = ("Could not automatically determine a suitable frequency " 1041 "at which to calculate the sensitivity. The overall " 1042 "sensitivity will not be recalculated.") 1043 raise ValueError(msg) 1044 1045 freq, gain = self._get_overall_sensitivity_and_gain( 1046 output=unit, frequency=float(frequency)) 1047 1048 self.instrument_sensitivity.value = gain 1049 self.instrument_sensitivity.frequency = freq 1050 1051 def _get_overall_sensitivity_and_gain( 1052 self, frequency=None, output='VEL'): 1053 """ 1054 Get the overall sensitivity and gain from stages 1 to N. 1055 1056 Returns the overall sensitivity frequency and gain, which can be 1057 used to create stage 0. 1058 1059 :type output: str 1060 :param output: Output units. One of: 1061 1062 ``"DISP"`` 1063 displacement, output unit is meters 1064 ``"VEL"`` 1065 velocity, output unit is meters/second 1066 ``"ACC"`` 1067 acceleration, output unit is meters/second**2 1068 1069 :type frequency: float 1070 :param frequency: Frequency to calculate overall sensitivity for in 1071 Hertz. Defaults to normalization frequency of stage 1. 1072 :rtype: :tuple: ( float, float ) 1073 :returns: frequency and gain at frequency. 1074 """ 1075 if frequency is None: 1076 # XXX is this safe enough, or should we lookup the stage sequence 1077 # XXX number explicitly? 1078 frequency = self.response_stages[0].normalization_frequency 1079 response_at_frequency = self._call_eval_resp_for_frequencies( 1080 frequencies=[frequency], output=output, 1081 hide_sensitivity_mismatch_warning=True)[0][0] 1082 overall_sensitivity = abs(response_at_frequency) 1083 return frequency, overall_sensitivity 1084 1085 def _call_eval_resp_for_frequencies( 1086 self, frequencies, output="VEL", start_stage=None, 1087 end_stage=None, hide_sensitivity_mismatch_warning=False): 1088 """ 1089 Returns frequency response for given frequencies using evalresp. 1090 1091 Also returns the overall sensitivity frequency and its gain. 1092 1093 :type frequencies: list of float 1094 :param frequencies: Discrete frequencies to calculate response for. 1095 :type output: str 1096 :param output: Output units. One of: 1097 1098 ``"DISP"`` 1099 displacement, output unit is meters 1100 ``"VEL"`` 1101 velocity, output unit is meters/second 1102 ``"ACC"`` 1103 acceleration, output unit is meters/second**2 1104 1105 :type start_stage: int, optional 1106 :param start_stage: Stage sequence number of first stage that will be 1107 used (disregarding all earlier stages). 1108 :type end_stage: int, optional 1109 :param end_stage: Stage sequence number of last stage that will be 1110 used (disregarding all later stages). 1111 :type hide_sensitivity_mismatch_warning: bool 1112 :param hide_sensitivity_mismatch_warning: Hide the evalresp warning 1113 that computed and reported sensitivities don't match. 1114 :rtype: :tuple: ( :class:`numpy.ndarray`, chan ) 1115 :returns: frequency response at requested frequencies 1116 """ 1117 if not self.response_stages: 1118 msg = ("Can not use evalresp on response with no response " 1119 "stages.") 1120 raise ObsPyException(msg) 1121 1122 import obspy.signal.evrespwrapper as ew 1123 from obspy.signal.headers import clibevresp 1124 1125 out_units = output.upper() 1126 if out_units not in ("DISP", "VEL", "ACC"): 1127 msg = ("requested output is '%s' but must be one of 'DISP', 'VEL' " 1128 "or 'ACC'") % output 1129 raise ValueError(msg) 1130 1131 frequencies = np.asarray(frequencies) 1132 1133 # Whacky. Evalresp uses a global variable and uses that to scale the 1134 # response if it encounters any unit that is not SI. 1135 scale_factor = [1.0] 1136 1137 def get_unit_mapping(key): 1138 try: 1139 key = key.upper() 1140 except Exception: 1141 pass 1142 units_mapping = { 1143 "M": ew.ENUM_UNITS["DIS"], 1144 "NM": ew.ENUM_UNITS["DIS"], 1145 "CM": ew.ENUM_UNITS["DIS"], 1146 "MM": ew.ENUM_UNITS["DIS"], 1147 "M/S": ew.ENUM_UNITS["VEL"], 1148 "M/SEC": ew.ENUM_UNITS["VEL"], 1149 "NM/S": ew.ENUM_UNITS["VEL"], 1150 "NM/SEC": ew.ENUM_UNITS["VEL"], 1151 "CM/S": ew.ENUM_UNITS["VEL"], 1152 "CM/SEC": ew.ENUM_UNITS["VEL"], 1153 "MM/S": ew.ENUM_UNITS["VEL"], 1154 "MM/SEC": ew.ENUM_UNITS["VEL"], 1155 "M/S**2": ew.ENUM_UNITS["ACC"], 1156 "M/(S**2)": ew.ENUM_UNITS["ACC"], 1157 "M/SEC**2": ew.ENUM_UNITS["ACC"], 1158 "M/(SEC**2)": ew.ENUM_UNITS["ACC"], 1159 "M/S/S": ew.ENUM_UNITS["ACC"], 1160 "NM/S**2": ew.ENUM_UNITS["ACC"], 1161 "NM/(S**2)": ew.ENUM_UNITS["ACC"], 1162 "NM/SEC**2": ew.ENUM_UNITS["ACC"], 1163 "NM/(SEC**2)": ew.ENUM_UNITS["ACC"], 1164 "CM/S**2": ew.ENUM_UNITS["ACC"], 1165 "CM/(S**2)": ew.ENUM_UNITS["ACC"], 1166 "CM/SEC**2": ew.ENUM_UNITS["ACC"], 1167 "CM/(SEC**2)": ew.ENUM_UNITS["ACC"], 1168 "MM/S**2": ew.ENUM_UNITS["ACC"], 1169 "MM/(S**2)": ew.ENUM_UNITS["ACC"], 1170 "MM/SEC**2": ew.ENUM_UNITS["ACC"], 1171 "MM/(SEC**2)": ew.ENUM_UNITS["ACC"], 1172 # Evalresp internally treats strain as displacement. 1173 "M/M": ew.ENUM_UNITS["DIS"], 1174 "M**3/M**3": ew.ENUM_UNITS["DIS"], 1175 "V": ew.ENUM_UNITS["VOLTS"], 1176 "VOLT": ew.ENUM_UNITS["VOLTS"], 1177 "VOLTS": ew.ENUM_UNITS["VOLTS"], 1178 # This is weird, but evalresp appears to do the same. 1179 "V/M": ew.ENUM_UNITS["VOLTS"], 1180 "COUNT": ew.ENUM_UNITS["COUNTS"], 1181 "COUNTS": ew.ENUM_UNITS["COUNTS"], 1182 "T": ew.ENUM_UNITS["TESLA"], 1183 "PA": ew.ENUM_UNITS["PRESSURE"], 1184 "PASCAL": ew.ENUM_UNITS["PRESSURE"], 1185 "PASCALS": ew.ENUM_UNITS["PRESSURE"], 1186 "MBAR": ew.ENUM_UNITS["PRESSURE"]} 1187 if key not in units_mapping: 1188 if key is not None: 1189 msg = ("The unit '%s' is not known to ObsPy. It will be " 1190 "assumed to be displacement for the calculations. " 1191 "This mostly does the right thing but please " 1192 "proceed with caution.") % key 1193 warnings.warn(msg) 1194 value = ew.ENUM_UNITS["DIS"] 1195 else: 1196 value = units_mapping[key] 1197 1198 # Scale factor with the same logic as evalresp. 1199 if key in ["CM/S**2", "CM/S", "CM/SEC", "CM"]: 1200 scale_factor[0] = 1.0E2 1201 elif key in ["MM/S**2", "MM/S", "MM/SEC", "MM"]: 1202 scale_factor[0] = 1.0E3 1203 elif key in ["NM/S**2", "NM/S", "NM/SEC", "NM"]: 1204 scale_factor[0] = 1.0E9 1205 1206 return value 1207 1208 all_stages = defaultdict(list) 1209 1210 for stage in self.response_stages: 1211 # optionally select only stages as requested by user 1212 if start_stage is not None: 1213 if stage.stage_sequence_number < start_stage: 1214 continue 1215 if end_stage is not None: 1216 if stage.stage_sequence_number > end_stage: 1217 continue 1218 all_stages[stage.stage_sequence_number].append(stage) 1219 1220 stage_lengths = set(map(len, all_stages.values())) 1221 if len(stage_lengths) != 1 or stage_lengths.pop() != 1: 1222 msg = "Each stage can only appear once." 1223 raise ValueError(msg) 1224 1225 stage_list = sorted(all_stages.keys()) 1226 1227 stage_objects = [] 1228 1229 # Attempt to fix some potentially faulty responses here. 1230 if 1 in all_stages and all_stages[1] and ( 1231 not all_stages[1][0].input_units or 1232 not all_stages[1][0].output_units): 1233 # Make a copy to not modify the original 1234 all_stages[1][0] = copy.deepcopy(all_stages[1][0]) 1235 # Some stages 1 are just the sensitivity and as thus don't store 1236 # input and output units in for example StationXML. In these cases 1237 # try to guess it from the overall sensitivity or stage 2. 1238 if not all_stages[1][0].input_units: 1239 if self.instrument_sensitivity.input_units: 1240 all_stages[1][0].input_units = \ 1241 self.instrument_sensitivity.input_units 1242 msg = "Set the input units of stage 1 to the overall " \ 1243 "input units." 1244 warnings.warn(msg) 1245 if not all_stages[1][0].output_units: 1246 if max(all_stages.keys()) == 1 and \ 1247 self.instrument_sensitivity.output_units: 1248 all_stages[1][0].output_units = \ 1249 self.instrument_sensitivity.output_units 1250 msg = "Set the output units of stage 1 to the overall " \ 1251 "output units." 1252 warnings.warn(msg) 1253 if 2 in all_stages and all_stages[2] and \ 1254 all_stages[2][0].input_units: 1255 all_stages[1][0].output_units = \ 1256 all_stages[2][0].input_units 1257 msg = "Set the output units of stage 1 to the input " \ 1258 "units of stage 2." 1259 warnings.warn(msg) 1260 1261 for stage_number in stage_list: 1262 st = ew.Stage() 1263 st.sequence_no = stage_number 1264 1265 stage_blkts = [] 1266 1267 blockette = all_stages[stage_number][0] 1268 1269 # Write the input and output units. 1270 st.input_units = get_unit_mapping(blockette.input_units) 1271 st.output_units = get_unit_mapping(blockette.output_units) 1272 1273 if isinstance(blockette, PolesZerosResponseStage): 1274 blkt = ew.Blkt() 1275 # Map the transfer function type. 1276 transfer_fct_mapping = { 1277 "LAPLACE (RADIANS/SECOND)": "LAPLACE_PZ", 1278 "LAPLACE (HERTZ)": "ANALOG_PZ", 1279 "DIGITAL (Z-TRANSFORM)": "IIR_PZ"} 1280 blkt.type = ew.ENUM_FILT_TYPES[transfer_fct_mapping[ 1281 blockette.pz_transfer_function_type]] 1282 1283 # The blockette is a pole zero blockette. 1284 pz = blkt.blkt_info.pole_zero 1285 1286 pz.nzeros = len(blockette.zeros) 1287 pz.npoles = len(blockette.poles) 1288 pz.a0 = blockette.normalization_factor 1289 pz.a0_freq = blockette.normalization_frequency 1290 1291 # XXX: Find a better way to do this. 1292 poles = (ew.ComplexNumber * len(blockette.poles))() 1293 for i, value in enumerate(blockette.poles): 1294 poles[i].real = value.real 1295 poles[i].imag = value.imag 1296 1297 zeros = (ew.ComplexNumber * len(blockette.zeros))() 1298 for i, value in enumerate(blockette.zeros): 1299 zeros[i].real = value.real 1300 zeros[i].imag = value.imag 1301 1302 pz.poles = C.cast(C.pointer(poles), 1303 C.POINTER(ew.ComplexNumber)) 1304 pz.zeros = C.cast(C.pointer(zeros), 1305 C.POINTER(ew.ComplexNumber)) 1306 elif isinstance(blockette, CoefficientsTypeResponseStage): 1307 blkt = ew.Blkt() 1308 # This type can have either an FIR or an IIR response. If 1309 # the number of denominators is 0, it is a FIR. Otherwise 1310 # an IIR. 1311 1312 # FIR 1313 if len(blockette.denominator) == 0: 1314 if blockette.cf_transfer_function_type.lower() \ 1315 != "digital": 1316 msg = ("When no denominators are given it must " 1317 "be a digital FIR filter.") 1318 raise ValueError(msg) 1319 # Set the type to an asymmetric FIR blockette. 1320 blkt.type = ew.ENUM_FILT_TYPES["FIR_ASYM"] 1321 fir = blkt.blkt_info.fir 1322 fir.h0 = 1.0 1323 fir.ncoeffs = len(blockette.numerator) 1324 # XXX: Find a better way to do this. 1325 coeffs = (C.c_double * len(blockette.numerator))() 1326 for i, value in enumerate(blockette.numerator): 1327 coeffs[i] = float(value) 1328 fir.coeffs = C.cast(C.pointer(coeffs), 1329 C.POINTER(C.c_double)) 1330 # IIR 1331 else: 1332 blkt.type = ew.ENUM_FILT_TYPES["IIR_COEFFS"] 1333 coeff = blkt.blkt_info.coeff 1334 1335 coeff.h0 = 1.0 1336 coeff.nnumer = len(blockette.numerator) 1337 coeff.ndenom = len(blockette.denominator) 1338 1339 # XXX: Find a better way to do this. 1340 coeffs = (C.c_double * len(blockette.numerator))() 1341 for i, value in enumerate(blockette.numerator): 1342 coeffs[i] = float(value) 1343 coeff.numer = C.cast(C.pointer(coeffs), 1344 C.POINTER(C.c_double)) 1345 coeffs = (C.c_double * len(blockette.denominator))() 1346 for i, value in enumerate(blockette.denominator): 1347 coeffs[i] = float(value) 1348 coeff.denom = C.cast(C.pointer(coeffs), 1349 C.POINTER(C.c_double)) 1350 elif isinstance(blockette, ResponseListResponseStage): 1351 blkt = ew.Blkt() 1352 blkt.type = ew.ENUM_FILT_TYPES["LIST"] 1353 1354 # Get values as numpy arrays. 1355 f = np.array([float(_i.frequency) 1356 for _i in blockette.response_list_elements], 1357 dtype=np.float64) 1358 amp = np.array([float(_i.amplitude) 1359 for _i in blockette.response_list_elements], 1360 dtype=np.float64) 1361 phase = np.array([ 1362 float(_i.phase) 1363 for _i in blockette.response_list_elements], 1364 dtype=np.float64) 1365 1366 # Sanity check. 1367 min_f = frequencies[frequencies > 0].min() 1368 max_f = frequencies.max() 1369 1370 min_f_avail = min(f) 1371 max_f_avail = max(f) 1372 1373 # Allow interpolation for at most two samples. 1374 _d = np.abs(np.diff(f)) 1375 _d = _d[_d > 0].min() * 2 1376 min_f_avail -= _d 1377 max_f_avail += _d 1378 1379 if min_f < min_f_avail or max_f > max_f_avail: 1380 msg = ( 1381 "Cannot calculate the response as it contains a " 1382 "response list stage with frequencies only from " 1383 "%.4f - %.4f Hz. You are requesting a response from " 1384 "%.4f - %.4f Hz.") 1385 raise ValueError(msg % (min_f_avail, max_f_avail, min_f, 1386 max_f)) 1387 1388 amp = scipy.interpolate.InterpolatedUnivariateSpline( 1389 f, amp, k=3)(frequencies) 1390 phase = scipy.interpolate.InterpolatedUnivariateSpline( 1391 f, phase, k=3)(frequencies) 1392 1393 # Set static offset to zero. 1394 amp[amp == 0] = 0 1395 phase[phase == 0] = 0 1396 1397 rl = blkt.blkt_info.list 1398 rl.nresp = len(frequencies) 1399 1400 _freq_c = (C.c_double * rl.nresp)() 1401 _amp_c = (C.c_double * rl.nresp)() 1402 _phase_c = (C.c_double * rl.nresp)() 1403 1404 _freq_c[:] = frequencies[:] 1405 _amp_c[:] = amp[:] 1406 _phase_c[:] = phase[:] 1407 1408 rl.freq = C.cast(C.pointer(_freq_c), 1409 C.POINTER(C.c_double)) 1410 rl.amp = C.cast(C.pointer(_amp_c), 1411 C.POINTER(C.c_double)) 1412 rl.phase = C.cast(C.pointer(_phase_c), 1413 C.POINTER(C.c_double)) 1414 1415 elif isinstance(blockette, FIRResponseStage): 1416 blkt = ew.Blkt() 1417 1418 if blockette.symmetry == "NONE": 1419 blkt.type = ew.ENUM_FILT_TYPES["FIR_ASYM"] 1420 if blockette.symmetry == "ODD": 1421 blkt.type = ew.ENUM_FILT_TYPES["FIR_SYM_1"] 1422 if blockette.symmetry == "EVEN": 1423 blkt.type = ew.ENUM_FILT_TYPES["FIR_SYM_2"] 1424 1425 # The blockette is a fir blockette 1426 fir = blkt.blkt_info.fir 1427 fir.h0 = 1.0 1428 fir.ncoeffs = len(blockette.coefficients) 1429 1430 # XXX: Find a better way to do this. 1431 coeffs = (C.c_double * len(blockette.coefficients))() 1432 for i, value in enumerate(blockette.coefficients): 1433 coeffs[i] = float(value) 1434 fir.coeffs = C.cast(C.pointer(coeffs), 1435 C.POINTER(C.c_double)) 1436 elif isinstance(blockette, PolynomialResponseStage): 1437 msg = ("PolynomialResponseStage not yet implemented. " 1438 "Please contact the developers.") 1439 raise NotImplementedError(msg) 1440 else: 1441 # Otherwise it could be a gain only stage. 1442 if blockette.stage_gain is not None and \ 1443 blockette.stage_gain_frequency is not None: 1444 blkt = None 1445 else: 1446 msg = "Type: %s." % str(type(blockette)) 1447 raise NotImplementedError(msg) 1448 1449 if blkt is not None: 1450 stage_blkts.append(blkt) 1451 1452 # Evalresp requires FIR and IIR blockettes to have decimation 1453 # values. Set the "unit decimation" values in case they are not 1454 # set. 1455 # 1456 # Only set it if there is a stage gain - otherwise evalresp 1457 # complains again. 1458 if isinstance(blockette, PolesZerosResponseStage) and \ 1459 blockette.stage_gain and \ 1460 None in set([ 1461 blockette.decimation_correction, 1462 blockette.decimation_delay, 1463 blockette.decimation_factor, 1464 blockette.decimation_input_sample_rate, 1465 blockette.decimation_offset]): 1466 # Don't modify the original object. 1467 blockette = copy.deepcopy(blockette) 1468 blockette.decimation_correction = 0.0 1469 blockette.decimation_delay = 0.0 1470 blockette.decimation_factor = 1 1471 blockette.decimation_offset = 0 1472 sr = self.get_sampling_rates() 1473 if sr and blockette.stage_sequence_number in sr and \ 1474 sr[blockette.stage_sequence_number][ 1475 "input_sampling_rate"]: 1476 blockette.decimation_input_sample_rate = \ 1477 self.get_sampling_rates()[ 1478 blockette.stage_sequence_number][ 1479 "input_sampling_rate"] 1480 # This branch get's large called for responses that only have a 1481 # a single stage. 1482 else: 1483 blockette.decimation_input_sample_rate = 1.0 1484 1485 # Parse the decimation if is given. 1486 decimation_values = set([ 1487 blockette.decimation_correction, 1488 blockette.decimation_delay, blockette.decimation_factor, 1489 blockette.decimation_input_sample_rate, 1490 blockette.decimation_offset]) 1491 if None in decimation_values: 1492 if len(decimation_values) != 1: 1493 msg = ("If a decimation is given, all values must " 1494 "be specified.") 1495 raise ValueError(msg) 1496 else: 1497 blkt = ew.Blkt() 1498 blkt.type = ew.ENUM_FILT_TYPES["DECIMATION"] 1499 decimation_blkt = blkt.blkt_info.decimation 1500 1501 # Evalresp does the same! 1502 if blockette.decimation_input_sample_rate == 0: 1503 decimation_blkt.sample_int = 0.0 1504 else: 1505 decimation_blkt.sample_int = \ 1506 1.0 / blockette.decimation_input_sample_rate 1507 1508 decimation_blkt.deci_fact = blockette.decimation_factor 1509 decimation_blkt.deci_offset = blockette.decimation_offset 1510 decimation_blkt.estim_delay = blockette.decimation_delay 1511 decimation_blkt.applied_corr = \ 1512 blockette.decimation_correction 1513 stage_blkts.append(blkt) 1514 1515 # Add the gain if it is available. 1516 if blockette.stage_gain is not None and \ 1517 blockette.stage_gain_frequency is not None: 1518 blkt = ew.Blkt() 1519 blkt.type = ew.ENUM_FILT_TYPES["GAIN"] 1520 gain_blkt = blkt.blkt_info.gain 1521 gain_blkt.gain = blockette.stage_gain 1522 gain_blkt.gain_freq = blockette.stage_gain_frequency 1523 stage_blkts.append(blkt) 1524 1525 if not stage_blkts: 1526 msg = "At least one blockette is needed for the stage." 1527 raise ValueError(msg) 1528 1529 # Attach the blockette chain to the stage. 1530 st.first_blkt = C.pointer(stage_blkts[0]) 1531 for _i in range(1, len(stage_blkts)): 1532 stage_blkts[_i - 1].next_blkt = C.pointer(stage_blkts[_i]) 1533 1534 stage_objects.append(st) 1535 1536 # Attach the instrument sensitivity as stage 0 at the end. 1537 st = ew.Stage() 1538 st.sequence_no = 0 1539 st.input_units = 0 1540 st.output_units = 0 1541 blkt = ew.Blkt() 1542 blkt.type = ew.ENUM_FILT_TYPES["GAIN"] 1543 gain_blkt = blkt.blkt_info.gain 1544 gain_blkt.gain = self.instrument_sensitivity.value 1545 # Set the sensitivity frequency - use 1.0 if not given. This is also 1546 # what evalresp does. 1547 gain_blkt.gain_freq = self.instrument_sensitivity.frequency \ 1548 if self.instrument_sensitivity.frequency else 0.0 1549 st.first_blkt = C.pointer(blkt) 1550 stage_objects.append(st) 1551 1552 chan = ew.Channel() 1553 if not stage_objects: 1554 msg = "At least one stage is needed." 1555 raise ValueError(msg) 1556 1557 # Attach the stage chain to the channel. 1558 chan.first_stage = C.pointer(stage_objects[0]) 1559 for _i in range(1, len(stage_objects)): 1560 stage_objects[_i - 1].next_stage = C.pointer(stage_objects[_i]) 1561 1562 chan.nstages = len(stage_objects) 1563 1564 # Evalresp will take care of setting it to the overall sensitivity. 1565 chan.sensit = 0.0 1566 chan.sensfreq = 0.0 1567 1568 output = np.empty(len(frequencies), dtype=np.complex128) 1569 out_units = C.c_char_p(out_units.encode('ascii', 'strict')) 1570 1571 # Set global variables 1572 if self.resource_id: 1573 clibevresp.curr_file.value = self.resource_id.encode('utf-8') 1574 else: 1575 clibevresp.curr_file.value = None 1576 1577 try: 1578 rc = clibevresp._obspy_check_channel(C.byref(chan)) 1579 if rc: 1580 e, m = ew.ENUM_ERROR_CODES[rc] 1581 raise e('check_channel: ' + m) 1582 rc = clibevresp._obspy_norm_resp( 1583 C.byref(chan), -1, 0, 1584 1 if hide_sensitivity_mismatch_warning else 0) 1585 if rc: 1586 e, m = ew.ENUM_ERROR_CODES[rc] 1587 raise e('norm_resp: ' + m) 1588 1589 rc = clibevresp._obspy_calc_resp(C.byref(chan), frequencies, 1590 len(frequencies), 1591 output, out_units, -1, 0, 0) 1592 if rc: 1593 e, m = ew.ENUM_ERROR_CODES[rc] 1594 raise e('calc_resp: ' + m) 1595 1596 # XXX: Check if this is really not needed. 1597 # output *= scale_factor[0] 1598 1599 finally: 1600 clibevresp.curr_file.value = None 1601 1602 return output, chan 1603 1604 def get_evalresp_response_for_frequencies( 1605 self, frequencies, output="VEL", start_stage=None, end_stage=None): 1606 """ 1607 Returns frequency response for given frequencies using evalresp. 1608 1609 :type frequencies: list of float 1610 :param frequencies: Discrete frequencies to calculate response for. 1611 :type output: str 1612 :param output: Output units. One of: 1613 1614 ``"DISP"`` 1615 displacement, output unit is meters 1616 ``"VEL"`` 1617 velocity, output unit is meters/second 1618 ``"ACC"`` 1619 acceleration, output unit is meters/second**2 1620 1621 :type start_stage: int, optional 1622 :param start_stage: Stage sequence number of first stage that will be 1623 used (disregarding all earlier stages). 1624 :type end_stage: int, optional 1625 :param end_stage: Stage sequence number of last stage that will be 1626 used (disregarding all later stages). 1627 :rtype: :class:`numpy.ndarray` 1628 :returns: frequency response at requested frequencies 1629 """ 1630 output, chan = self._call_eval_resp_for_frequencies( 1631 frequencies, output=output, start_stage=start_stage, 1632 end_stage=end_stage) 1633 return output 1634 1635 def get_evalresp_response(self, t_samp, nfft, output="VEL", 1636 start_stage=None, end_stage=None): 1637 """ 1638 Returns frequency response and corresponding frequencies using 1639 evalresp. 1640 1641 :type t_samp: float 1642 :param t_samp: time resolution (inverse frequency resolution) 1643 :type nfft: int 1644 :param nfft: Number of FFT points to use 1645 :type output: str 1646 :param output: Output units. One of: 1647 1648 ``"DISP"`` 1649 displacement, output unit is meters 1650 ``"VEL"`` 1651 velocity, output unit is meters/second 1652 ``"ACC"`` 1653 acceleration, output unit is meters/second**2 1654 1655 :type start_stage: int, optional 1656 :param start_stage: Stage sequence number of first stage that will be 1657 used (disregarding all earlier stages). 1658 :type end_stage: int, optional 1659 :param end_stage: Stage sequence number of last stage that will be 1660 used (disregarding all later stages). 1661 :rtype: tuple of two arrays 1662 :returns: frequency response and corresponding frequencies 1663 """ 1664 # Calculate the output frequencies. 1665 fy = 1 / (t_samp * 2.0) 1666 # start at zero to get zero for offset/ DC of fft 1667 # numpy 1.9 introduced a dtype kwarg 1668 try: 1669 freqs = np.linspace(0, fy, int(nfft // 2) + 1, dtype=np.float64) 1670 except Exception: 1671 freqs = np.linspace(0, fy, int(nfft // 2) + 1).astype(np.float64) 1672 1673 response = self.get_evalresp_response_for_frequencies( 1674 freqs, output=output, start_stage=start_stage, end_stage=end_stage) 1675 return response, freqs 1676 1677 def __str__(self): 1678 i_s = self.instrument_sensitivity 1679 if i_s: 1680 input_units = i_s.input_units \ 1681 if i_s.input_units else "UNKNOWN" 1682 input_units_description = i_s.input_units_description \ 1683 if i_s.input_units_description else "" 1684 output_units = i_s.output_units \ 1685 if i_s.output_units else "UNKNOWN" 1686 output_units_description = i_s.output_units_description \ 1687 if i_s.output_units_description else "" 1688 sensitivity = ("%g" % i_s.value) if i_s.value else "UNKNOWN" 1689 freq = ("%.3f" % i_s.frequency) if i_s.frequency else "UNKNOWN" 1690 else: 1691 input_units = "UNKNOWN" 1692 input_units_description = "" 1693 output_units = "UNKNOWN" 1694 output_units_description = "" 1695 sensitivity = "UNKNOWN" 1696 freq = "UNKNOWN" 1697 1698 ret = ( 1699 "Channel Response\n" 1700 "\tFrom {input_units} ({input_units_description}) to " 1701 "{output_units} ({output_units_description})\n" 1702 "\tOverall Sensitivity: {sensitivity} defined at {freq} Hz\n" 1703 "\t{stages} stages:\n{stage_desc}").format( 1704 input_units=input_units, 1705 input_units_description=input_units_description, 1706 output_units=output_units, 1707 output_units_description=output_units_description, 1708 sensitivity=sensitivity, 1709 freq=freq, 1710 stages=len(self.response_stages), 1711 stage_desc="\n".join( 1712 ["\t\tStage %i: %s from %s to %s," 1713 " gain: %s" % ( 1714 i.stage_sequence_number, i.__class__.__name__, 1715 i.input_units, i.output_units, 1716 ("%g" % i.stage_gain) if i.stage_gain else "UNKNOWN") 1717 for i in self.response_stages])) 1718 return ret 1719 1720 def _repr_pretty_(self, p, cycle): 1721 p.text(str(self)) 1722 1723 def plot(self, min_freq, output="VEL", start_stage=None, 1724 end_stage=None, label=None, axes=None, sampling_rate=None, 1725 unwrap_phase=False, plot_degrees=False, show=True, outfile=None): 1726 """ 1727 Show bode plot of instrument response. 1728 1729 :type min_freq: float 1730 :param min_freq: Lowest frequency to plot. 1731 :type output: str 1732 :param output: Output units. One of: 1733 1734 ``"DISP"`` 1735 displacement 1736 ``"VEL"`` 1737 velocity 1738 ``"ACC"`` 1739 acceleration 1740 1741 :type start_stage: int, optional 1742 :param start_stage: Stage sequence number of first stage that will be 1743 used (disregarding all earlier stages). 1744 :type end_stage: int, optional 1745 :param end_stage: Stage sequence number of last stage that will be 1746 used (disregarding all later stages). 1747 :type label: str 1748 :param label: Label string for legend. 1749 :type axes: list of 2 :class:`matplotlib.axes.Axes` 1750 :param axes: List/tuple of two axes instances to plot the 1751 amplitude/phase spectrum into. If not specified, a new figure is 1752 opened. 1753 :type sampling_rate: float 1754 :param sampling_rate: Manually specify sampling rate of time series. 1755 If not given it is attempted to determine it from the information 1756 in the individual response stages. Does not influence the spectra 1757 calculation, if it is not known, just provide the highest frequency 1758 that should be plotted times two. 1759 :type unwrap_phase: bool 1760 :param unwrap_phase: Set optional phase unwrapping using NumPy. 1761 :type plot_degrees: bool 1762 :param plot_degrees: if ``True`` plot bode in degrees 1763 :type show: bool 1764 :param show: Whether to show the figure after plotting or not. Can be 1765 used to do further customization of the plot before showing it. 1766 :type outfile: str 1767 :param outfile: Output file path to directly save the resulting image 1768 (e.g. ``"/tmp/image.png"``). Overrides the ``show`` option, image 1769 will not be displayed interactively. The given path/filename is 1770 also used to automatically determine the output format. Supported 1771 file formats depend on your matplotlib backend. Most backends 1772 support png, pdf, ps, eps and svg. Defaults to ``None``. 1773 1774 .. rubric:: Basic Usage 1775 1776 >>> from obspy import read_inventory 1777 >>> resp = read_inventory()[0][0][0].response 1778 >>> resp.plot(0.001, output="VEL") # doctest: +SKIP 1779 1780 .. plot:: 1781 1782 from obspy import read_inventory 1783 resp = read_inventory()[0][0][0].response 1784 resp.plot(0.001, output="VEL") 1785 """ 1786 import matplotlib.pyplot as plt 1787 from matplotlib.transforms import blended_transform_factory 1788 1789 # detect sampling rate from response stages 1790 if sampling_rate is None: 1791 for stage in self.response_stages[::-1]: 1792 if (stage.decimation_input_sample_rate is not None and 1793 stage.decimation_factor is not None): 1794 sampling_rate = (stage.decimation_input_sample_rate / 1795 stage.decimation_factor) 1796 break 1797 else: 1798 msg = ("Failed to autodetect sampling rate of channel from " 1799 "response stages. Please manually specify parameter " 1800 "`sampling_rate`") 1801 raise Exception(msg) 1802 if sampling_rate == 0: 1803 msg = "Can not plot response for channel with sampling rate `0`." 1804 raise ZeroSamplingRate(msg) 1805 1806 t_samp = 1.0 / sampling_rate 1807 nyquist = sampling_rate / 2.0 1808 nfft = int(sampling_rate / min_freq) 1809 1810 cpx_response, freq = self.get_evalresp_response( 1811 t_samp=t_samp, nfft=nfft, output=output, start_stage=start_stage, 1812 end_stage=end_stage) 1813 1814 if axes is not None: 1815 ax1, ax2 = axes 1816 fig = ax1.figure 1817 else: 1818 fig = plt.figure() 1819 ax1 = fig.add_subplot(211) 1820 ax2 = fig.add_subplot(212, sharex=ax1) 1821 1822 label_kwarg = {} 1823 if label is not None: 1824 label_kwarg['label'] = label 1825 1826 # plot amplitude response 1827 lw = 1.5 1828 lines = ax1.loglog(freq, abs(cpx_response), lw=lw, **label_kwarg) 1829 color = lines[0].get_color() 1830 if self.instrument_sensitivity: 1831 trans_above = blended_transform_factory(ax1.transData, 1832 ax1.transAxes) 1833 trans_right = blended_transform_factory(ax1.transAxes, 1834 ax1.transData) 1835 arrowprops = dict( 1836 arrowstyle="wedge,tail_width=1.4,shrink_factor=0.8", fc=color) 1837 bbox = dict(boxstyle="round", fc="w") 1838 if self.instrument_sensitivity.frequency: 1839 ax1.annotate("%.1g" % self.instrument_sensitivity.frequency, 1840 (self.instrument_sensitivity.frequency, 1.0), 1841 xytext=(self.instrument_sensitivity.frequency, 1842 1.1), 1843 xycoords=trans_above, textcoords=trans_above, 1844 ha="center", va="bottom", 1845 arrowprops=arrowprops, bbox=bbox) 1846 if self.instrument_sensitivity.value: 1847 ax1.annotate("%.1e" % self.instrument_sensitivity.value, 1848 (1.0, self.instrument_sensitivity.value), 1849 xytext=(1.05, self.instrument_sensitivity.value), 1850 xycoords=trans_right, textcoords=trans_right, 1851 ha="left", va="center", 1852 arrowprops=arrowprops, bbox=bbox) 1853 1854 # plot phase response 1855 phase = np.angle(cpx_response, deg=plot_degrees) 1856 if unwrap_phase and not plot_degrees: 1857 phase = np.unwrap(phase) 1858 ax2.semilogx(freq, phase, color=color, lw=lw) 1859 1860 # plot nyquist frequency 1861 for ax in (ax1, ax2): 1862 ax.axvline(nyquist, ls="--", color=color, lw=lw) 1863 1864 # only do adjustments if we initialized the figure in here 1865 if not axes: 1866 _adjust_bode_plot_figure(fig, show=False, 1867 plot_degrees=plot_degrees) 1868 1869 if outfile: 1870 fig.savefig(outfile) 1871 else: 1872 if show: 1873 plt.show() 1874 1875 return fig 1876 1877 def get_paz(self): 1878 """ 1879 Get Poles and Zeros stage. 1880 1881 Prints a warning if more than one poles and zeros stage is found. 1882 Raises if no poles and zeros stage is found. 1883 1884 :rtype: :class:`PolesZerosResponseStage` 1885 :returns: Poles and Zeros response stage. 1886 """ 1887 paz = [deepcopy(stage) for stage in self.response_stages 1888 if isinstance(stage, PolesZerosResponseStage)] 1889 if len(paz) == 0: 1890 msg = "No PolesZerosResponseStage found." 1891 raise Exception(msg) 1892 elif len(paz) > 1: 1893 msg = ("More than one PolesZerosResponseStage encountered. " 1894 "Returning first one found.") 1895 warnings.warn(msg) 1896 return paz[0] 1897 1898 def get_sacpz(self): 1899 """ 1900 Returns SACPZ ASCII text representation of Response. 1901 1902 :rtype: str 1903 :returns: Textual SACPZ representation of response. 1904 """ 1905 # extract paz 1906 paz = self.get_paz() 1907 return paz_to_sacpz_string(paz, self.instrument_sensitivity) 1908 1909 @classmethod 1910 def from_paz(cls, zeros, poles, stage_gain, 1911 stage_gain_frequency=1.0, input_units='M/S', 1912 output_units='VOLTS', normalization_frequency=1.0, 1913 pz_transfer_function_type='LAPLACE (RADIANS/SECOND)', 1914 normalization_factor=1.0): 1915 """ 1916 Convert poles and zeros lists into a single-stage response. 1917 1918 Takes in lists of complex poles and zeros and returns a Response with 1919 those values defining its only stage. Most of the optional parameters 1920 defined here are from 1921 :class:`~obspy.core.inventory.response.PolesZerosResponseStage`. 1922 1923 :type zeros: list of complex 1924 :param zeros: All zeros of the response to be defined. 1925 :type poles: list of complex 1926 :param poles: All poles of the response to be defined. 1927 :type stage_gain: float 1928 :param stage_gain: The gain value of the response [sensitivity] 1929 :returns: new Response instance with given P-Z values 1930 """ 1931 pzstage = PolesZerosResponseStage( 1932 stage_sequence_number=1, stage_gain=stage_gain, 1933 stage_gain_frequency=stage_gain_frequency, 1934 input_units=input_units, output_units=output_units, 1935 pz_transfer_function_type=pz_transfer_function_type, 1936 normalization_frequency=normalization_frequency, 1937 zeros=zeros, poles=poles, 1938 normalization_factor=normalization_factor) 1939 sens = InstrumentSensitivity(value=stage_gain, 1940 frequency=stage_gain_frequency, 1941 input_units=input_units, 1942 output_units=output_units) 1943 resp = cls(instrument_sensitivity=sens, response_stages=[pzstage]) 1944 resp.recalculate_overall_sensitivity() 1945 return resp 1946 1947 1948def paz_to_sacpz_string(paz, instrument_sensitivity): 1949 """ 1950 Returns SACPZ ASCII text representation of Response. 1951 1952 :type paz: :class:`PolesZerosResponseStage` 1953 :param paz: Poles and Zeros response information 1954 :type instrument_sensitivity: :class:`InstrumentSensitivity` 1955 :param paz: Overall instrument sensitivity of response 1956 :rtype: str 1957 :returns: Textual SACPZ representation of poles and zeros response stage. 1958 """ 1959 # assemble output string 1960 out = [] 1961 out.append("ZEROS %i" % len(paz.zeros)) 1962 for c in paz.zeros: 1963 out.append(" %+.6e %+.6e" % (c.real, c.imag)) 1964 out.append("POLES %i" % len(paz.poles)) 1965 for c in paz.poles: 1966 out.append(" %+.6e %+.6e" % (c.real, c.imag)) 1967 constant = paz.normalization_factor * instrument_sensitivity.value 1968 out.append("CONSTANT %.6e" % constant) 1969 return "\n".join(out) 1970 1971 1972class InstrumentSensitivity(ComparingObject): 1973 """ 1974 From the StationXML Definition: 1975 The total sensitivity for a channel, representing the complete 1976 acquisition system expressed as a scalar. Equivalent to SEED stage 0 1977 gain with (blockette 58) with the ability to specify a frequency range. 1978 1979 Sensitivity and frequency ranges. The FrequencyRangeGroup is an optional 1980 construct that defines a pass band in Hertz (FrequencyStart and 1981 FrequencyEnd) in which the SensitivityValue is valid within the number of 1982 decibels specified in FrequencyDBVariation. 1983 """ 1984 def __init__(self, value, frequency, input_units, 1985 output_units, input_units_description=None, 1986 output_units_description=None, frequency_range_start=None, 1987 frequency_range_end=None, frequency_range_db_variation=None): 1988 """ 1989 :type value: float 1990 :param value: Complex type for sensitivity and frequency ranges. 1991 This complex type can be used to represent both overall 1992 sensitivities and individual stage gains. The FrequencyRangeGroup 1993 is an optional construct that defines a pass band in Hertz ( 1994 FrequencyStart and FrequencyEnd) in which the SensitivityValue is 1995 valid within the number of decibels specified in 1996 FrequencyDBVariation. 1997 :type frequency: float 1998 :param frequency: Complex type for sensitivity and frequency 1999 ranges. This complex type can be used to represent both overall 2000 sensitivities and individual stage gains. The FrequencyRangeGroup 2001 is an optional construct that defines a pass band in Hertz ( 2002 FrequencyStart and FrequencyEnd) in which the SensitivityValue is 2003 valid within the number of decibels specified in 2004 FrequencyDBVariation. 2005 :param input_units: string 2006 :param input_units: The units of the data as input from the 2007 perspective of data acquisition. After correcting data for this 2008 response, these would be the resulting units. 2009 Name of units, e.g. "M/S", "V", "PA". 2010 :param input_units_description: string, optional 2011 :param input_units_description: The units of the data as input from the 2012 perspective of data acquisition. After correcting data for this 2013 response, these would be the resulting units. 2014 Description of units, e.g. "Velocity in meters per second", 2015 "Volts", "Pascals". 2016 :param output_units: string 2017 :param output_units: The units of the data as output from the 2018 perspective of data acquisition. These would be the units of the 2019 data prior to correcting for this response. 2020 Name of units, e.g. "M/S", "V", "PA". 2021 :type output_units_description: str, optional 2022 :param output_units_description: The units of the data as output from 2023 the perspective of data acquisition. These would be the units of 2024 the data prior to correcting for this response. 2025 Description of units, e.g. "Velocity in meters per second", 2026 "Volts", "Pascals". 2027 :type frequency_range_start: float, optional 2028 :param frequency_range_start: Start of the frequency range for which 2029 the SensitivityValue is valid within the dB variation specified. 2030 :type frequency_range_end: float, optional 2031 :param frequency_range_end: End of the frequency range for which the 2032 SensitivityValue is valid within the dB variation specified. 2033 :type frequency_range_db_variation: float, optional 2034 :param frequency_range_db_variation: Variation in decibels within the 2035 specified range. 2036 """ 2037 self.value = value 2038 self.frequency = frequency 2039 self.input_units = input_units 2040 self.input_units_description = input_units_description 2041 self.output_units = output_units 2042 self.output_units_description = output_units_description 2043 self.frequency_range_start = frequency_range_start 2044 self.frequency_range_end = frequency_range_end 2045 self.frequency_range_db_variation = frequency_range_db_variation 2046 2047 def __str__(self): 2048 ret = ("Instrument Sensitivity:\n" 2049 "\tValue: {value}\n" 2050 "\tFrequency: {frequency}\n" 2051 "\tInput units: {input_units}\n" 2052 "\tInput units description: {input_units_description}\n" 2053 "\tOutput units: {output_units}\n" 2054 "\tOutput units description: {output_units_description}\n") 2055 ret = ret.format(**self.__dict__) 2056 return ret 2057 2058 def _repr_pretty_(self, p, cycle): 2059 p.text(str(self)) 2060 2061 2062# XXX duplicated code, PolynomialResponseStage could probably be implemented by 2063# XXX inheriting from both InstrumentPolynomial and ResponseStage 2064class InstrumentPolynomial(ComparingObject): 2065 """ 2066 From the StationXML Definition: 2067 The total sensitivity for a channel, representing the complete 2068 acquisition system expressed as a polynomial. Equivalent to SEED stage 2069 0 polynomial (blockette 62). 2070 """ 2071 def __init__(self, input_units, output_units, 2072 frequency_lower_bound, 2073 frequency_upper_bound, approximation_lower_bound, 2074 approximation_upper_bound, maximum_error, coefficients, 2075 approximation_type='MACLAURIN', resource_id=None, name=None, 2076 input_units_description=None, 2077 output_units_description=None, description=None): 2078 """ 2079 :type approximation_type: str 2080 :param approximation_type: Approximation type. Currently restricted to 2081 'MACLAURIN' by StationXML definition. 2082 :type frequency_lower_bound: float 2083 :param frequency_lower_bound: Lower frequency bound. 2084 :type frequency_upper_bound: float 2085 :param frequency_upper_bound: Upper frequency bound. 2086 :type approximation_lower_bound: float 2087 :param approximation_lower_bound: Lower bound of approximation. 2088 :type approximation_upper_bound: float 2089 :param approximation_upper_bound: Upper bound of approximation. 2090 :type maximum_error: float 2091 :param maximum_error: Maximum error. 2092 :type coefficients: list of floats 2093 :param coefficients: List of polynomial coefficients. 2094 :param input_units: string 2095 :param input_units: The units of the data as input from the 2096 perspective of data acquisition. After correcting data for this 2097 response, these would be the resulting units. 2098 Name of units, e.g. "M/S", "V", "PA". 2099 :param output_units: string 2100 :param output_units: The units of the data as output from the 2101 perspective of data acquisition. These would be the units of the 2102 data prior to correcting for this response. 2103 Name of units, e.g. "M/S", "V", "PA". 2104 :type resource_id: str 2105 :param resource_id: This field contains a string that should serve as a 2106 unique resource identifier. This identifier can be interpreted 2107 differently depending on the data center/software that generated 2108 the document. Also, we recommend to use something like 2109 GENERATOR:Meaningful ID. As a common behavior equipment with the 2110 same ID should contains the same information/be derived from the 2111 same base instruments. 2112 :type name: str 2113 :param name: A name given to the filter stage. 2114 :param input_units_description: string, optional 2115 :param input_units_description: The units of the data as input from the 2116 perspective of data acquisition. After correcting data for this 2117 response, these would be the resulting units. 2118 Description of units, e.g. "Velocity in meters per second", 2119 "Volts", "Pascals". 2120 :type output_units_description: str, optional 2121 :param output_units_description: The units of the data as output from 2122 the perspective of data acquisition. These would be the units of 2123 the data prior to correcting for this response. 2124 Description of units, e.g. "Velocity in meters per second", 2125 "Volts", "Pascals". 2126 :type description: str, optional 2127 :param description: A short description of of the filter. 2128 """ 2129 self.input_units = input_units 2130 self.output_units = output_units 2131 self.input_units_description = input_units_description 2132 self.output_units_description = output_units_description 2133 self.resource_id = resource_id 2134 self.name = name 2135 self.description = description 2136 self._approximation_type = approximation_type 2137 self.frequency_lower_bound = frequency_lower_bound 2138 self.frequency_upper_bound = frequency_upper_bound 2139 self.approximation_lower_bound = approximation_lower_bound 2140 self.approximation_upper_bound = approximation_upper_bound 2141 self.maximum_error = maximum_error 2142 self.coefficients = coefficients 2143 2144 @property 2145 def approximation_type(self): 2146 return self._approximation_type 2147 2148 @approximation_type.setter 2149 def approximation_type(self, value): 2150 value = str(value).upper() 2151 allowed = ("MACLAURIN",) 2152 if value not in allowed: 2153 msg = ("Value '%s' for polynomial response approximation type not " 2154 "allowed. Possible values are: '%s'") 2155 msg = msg % (value, "', '".join(allowed)) 2156 raise ValueError(msg) 2157 self._approximation_type = value 2158 2159 2160class FilterCoefficient(FloatWithUncertainties): 2161 """ 2162 A filter coefficient. 2163 """ 2164 def __init__(self, value, number=None): 2165 """ 2166 :type value: float 2167 :param value: The actual value of the coefficient 2168 :type number: int, optional 2169 :param number: Number to indicate the position of the coefficient. 2170 """ 2171 super(FilterCoefficient, self).__init__(value) 2172 self.number = number 2173 2174 @property 2175 def number(self): 2176 return self._number 2177 2178 @number.setter 2179 def number(self, value): 2180 if value is not None: 2181 value = int(value) 2182 self._number = value 2183 2184 2185class CoefficientWithUncertainties(FloatWithUncertainties): 2186 """ 2187 A coefficient with optional uncertainties. 2188 """ 2189 def __init__(self, value, number=None, lower_uncertainty=None, 2190 upper_uncertainty=None): 2191 """ 2192 :type value: float 2193 :param value: The actual value of the coefficient 2194 :type number: int, optional 2195 :param number: Number to indicate the position of the coefficient. 2196 :type lower_uncertainty: float 2197 :param lower_uncertainty: Lower uncertainty (aka minusError) 2198 :type upper_uncertainty: float 2199 :param upper_uncertainty: Upper uncertainty (aka plusError) 2200 """ 2201 super(CoefficientWithUncertainties, self).__init__( 2202 value, lower_uncertainty=lower_uncertainty, 2203 upper_uncertainty=upper_uncertainty) 2204 self.number = number 2205 2206 @property 2207 def number(self): 2208 return self._number 2209 2210 @number.setter 2211 def number(self, value): 2212 if value is not None: 2213 value = int(value) 2214 self._number = value 2215 2216 2217def _adjust_bode_plot_figure(fig, plot_degrees=False, grid=True, show=True): 2218 """ 2219 Helper function to do final adjustments to Bode plot figure. 2220 """ 2221 import matplotlib.pyplot as plt 2222 # make more room in between subplots for the ylabel of right plot 2223 fig.subplots_adjust(hspace=0.02, top=0.87, right=0.82) 2224 ax1, ax2 = fig.axes[:2] 2225 # workaround for older matplotlib versions 2226 try: 2227 ax1.legend(loc="lower center", ncol=3, fontsize='small') 2228 except TypeError: 2229 leg_ = ax1.legend(loc="lower center", ncol=3) 2230 leg_.prop.set_size("small") 2231 plt.setp(ax1.get_xticklabels(), visible=False) 2232 ax1.set_ylabel('Amplitude') 2233 minmax1 = ax1.get_ylim() 2234 ax1.set_ylim(top=minmax1[1] * 5) 2235 ax1.grid(True) 2236 ax2.set_xlabel('Frequency [Hz]') 2237 if plot_degrees: 2238 # degrees bode plot 2239 ax2.set_ylabel('Phase [degrees]') 2240 ax2.set_yticks(np.arange(-180, 180, 30)) 2241 ax2.set_yticklabels(np.arange(-180, 180, 30)) 2242 ax2.set_ylim(-180, 180) 2243 else: 2244 # radian bode plot 2245 plt.setp(ax2.get_yticklabels()[-1], visible=False) 2246 ax2.set_ylabel('Phase [rad]') 2247 minmax2 = ax2.yaxis.get_data_interval() 2248 yticks2 = np.arange(minmax2[0] - minmax2[0] % (pi / 2), 2249 minmax2[1] - minmax2[1] % (pi / 2) + pi, pi / 2) 2250 ax2.set_yticks(yticks2) 2251 ax2.set_yticklabels([_pitick2latex(x) for x in yticks2]) 2252 ax2.grid(True) 2253 2254 if show: 2255 plt.show() 2256 2257 2258def _pitick2latex(x): 2259 """ 2260 Helper function to convert a float that is a multiple of pi/2 2261 to a latex string. 2262 """ 2263 # safety check, if no multiple of pi/2 return normal representation 2264 if x % (pi / 2) != 0: 2265 return "%#.3g" % x 2266 string = "$" 2267 if x < 0: 2268 string += "-" 2269 if x / pi % 1 == 0: 2270 x = abs(int(x / pi)) 2271 if x == 0: 2272 return "$0$" 2273 elif x == 1: 2274 x = "" 2275 string += r"%s\pi$" % x 2276 else: 2277 x = abs(int(2 * x / pi)) 2278 if x == 1: 2279 x = "" 2280 string += r"\frac{%s\pi}{2}$" % x 2281 return string 2282 2283 2284if __name__ == '__main__': 2285 import doctest 2286 doctest.testmod(exclude_empty=True) 2287
