# -*- coding: utf-8 -*-
# ------------------------------------------------------------------
# Filename: calibration.py
# Purpose: Functions for relative calibration (e.g. Huddle test calibration)
# Author: Felix Bernauer, Simon Kremers
# Email: bernauer@geophysik.uni-muenchen.de
#
# Copyright (C) 2011 Felix Bernauer, Simon Kremers
# --------------------------------------------------------------------
"""
Functions for relative calibration.
:copyright:
The ObsPy Development Team (devs@obspy.org)
:license:
GNU Lesser General Public License, Version 3
(https://www.gnu.org/copyleft/lesser.html)
"""
import numpy as np
from obspy.core.stream import Stream
from obspy.core.trace import Trace
from obspy.io.gse2.paz import read_paz
from obspy.signal.invsim import paz_to_freq_resp
from obspy.signal.konnoohmachismoothing import konno_ohmachi_smoothing
from obspy.signal.util import next_pow_2
[docs]def rel_calib_stack(st1, st2, calib_file, window_len, overlap_frac=0.5,
smooth=0, save_data=True):
"""
Method for relative calibration of sensors using a sensor with known
transfer function
:param st1: Stream or Trace object, (known)
:param st2: Stream or Trace object, (unknown)
:type calib_file: str or dict
:param calib_file: file name of calibration file containing the PAZ of the
known instrument in GSE2 standard or a dictionary with poles and zeros
information (with keys ``'poles'``, ``'zeros'`` and ``'sensitivity'``).
:type window_len: float
:param window_len: length of sliding window in seconds
:type overlap_frac: float
:param overlap_frac: fraction of overlap, defaults to fifty percent (0.5)
:type smooth: float
:param smooth: variable that defines if the Konno-Ohmachi taper is used or
not. default = 0 -> no taper generally used in geopsy: smooth = 40
:type save_data: bool
:param save_data: Whether or not to save the result to a file. If True, two
output files will be created:
* The new response in station_name.window_length.resp
* The ref response in station_name.refResp
Defaults to True
:returns: frequency, amplitude and phase spectrum
implemented after rel_calib_stack.c by M.Ohrnberger and J.Wassermann.
"""
# transform given trace objects to streams
if isinstance(st1, Trace):
st1 = Stream([st1])
if isinstance(st2, Trace):
st2 = Stream([st2])
# check if sampling rate and trace length is the same
if st1[0].stats.npts != st2[0].stats.npts:
msg = "Traces don't have the same length!"
raise ValueError(msg)
elif st1[0].stats.sampling_rate != st2[0].stats.sampling_rate:
msg = "Traces don't have the same sampling rate!"
raise ValueError(msg)
else:
sampfreq = st1[0].stats.sampling_rate
# read waveforms
tr1 = st1[0].data.astype(np.float64)
tr2 = st2[0].data.astype(np.float64)
# get window length, nfft and frequency step
ndat = int(window_len * sampfreq)
nfft = next_pow_2(ndat)
# read calib file and calculate response function
gg, _freq = _calc_resp(calib_file, nfft, sampfreq)
# calculate number of windows and overlap
noverlap = nfft * overlap_frac
auto, _freq, _t = \
spectral_helper(tr1, tr1, NFFT=nfft, Fs=sampfreq, noverlap=noverlap)
cross, freq, _t = \
spectral_helper(tr2, tr1, NFFT=nfft, Fs=sampfreq, noverlap=noverlap)
# get number of windows that were actually computed inside FFT routine
nwin = auto.shape[1]
res = (cross / auto).sum(axis=1) * gg
# The first item might be zero. Problems with phase calculations.
res = res[1:]
freq = freq[1:]
gg = gg[1:]
res /= nwin
# apply Konno-Ohmachi smoothing taper if chosen
if smooth > 0:
# Write in one matrix for performance reasons.
spectra = np.empty((2, len(res.real)))
spectra[0] = res.real
spectra[1] = res.imag
new_spectra = \
konno_ohmachi_smoothing(spectra, freq, bandwidth=smooth, count=1,
max_memory_usage=1024, normalize=True)
res.real = new_spectra[0]
res.imag = new_spectra[1]
amp = np.abs(res)
# include phase unwrapping
phase = np.unwrap(np.angle(res)) # + 2.0 * np.pi
ra = np.abs(gg)
rpha = np.unwrap(np.angle(gg))
if save_data:
trans_new = (st2[0].stats.station + "." + st2[0].stats.channel +
"." + str(window_len) + ".resp")
trans_ref = st1[0].stats.station + ".refResp"
# Create empty array for easy saving
temp = np.empty((len(freq), 3))
temp[:, 0] = freq
temp[:, 1] = amp
temp[:, 2] = phase
np.savetxt(trans_new, temp, fmt='%.10f')
temp[:, 1] = ra
temp[:, 2] = rpha
np.savetxt(trans_ref, temp, fmt='%.10f')
return freq, amp, phase
[docs]def _calc_resp(calfile, nfft, sampfreq):
"""
Calculate transfer function of known system.
:type calfile: str
:param calfile: file containing poles, zeros and scale factor for known
system or a dictionary with poles and zeros information (with keys
``'poles'``, ``'zeros'`` and ``'sensitivity'``).
:returns: complex transfer function, array of frequencies
"""
# test if calfile is a paz dict
if isinstance(calfile, dict):
paz = calfile
# or read paz file if a filename is specified
else:
paz = dict()
paz['poles'], paz['zeros'], paz['sensitivity'] = read_paz(calfile)
# calculate transfer function
h, f = paz_to_freq_resp(paz['poles'], paz['zeros'], paz['sensitivity'],
1.0 / sampfreq, nfft, freq=True)
return h, f
# A modified copy of the Matplotlib 0.99.1.1 method spectral_helper found in
# .../matlab/mlab.py.
# Some function were changed to avoid additional dependencies. Included here as
# it is essential for the above rel_calib_stack function and only present in
# recent matplotlib versions.
# This is a helper function that implements the commonality between the
# psd, csd, and spectrogram. It is *NOT* meant to be used outside of mlab
[docs]def spectral_helper(x, y, NFFT=256, Fs=2, noverlap=0, pad_to=None, # noqa
sides='default', scale_by_freq=None):
# The checks for if y is x are so that we can use the same function to
# implement the core of psd(), csd(), and spectrogram() without doing
# extra calculations. We return the unaveraged Pxy, freqs, and t.
same_data = y is x
# Make sure we're dealing with a NumPy array. If y and x were the same
# object to start with, keep them that way
x = np.asarray(x)
if not same_data:
y = np.asarray(y)
# zero pad x and y up to NFFT if they are shorter than NFFT
if len(x) < NFFT:
n = len(x)
x = np.resize(x, (NFFT,))
x[n:] = 0
if not same_data and len(y) < NFFT:
n = len(y)
y = np.resize(y, (NFFT,))
y[n:] = 0
if pad_to is None:
pad_to = NFFT
if scale_by_freq is None:
scale_by_freq = True
# For real x, ignore the negative frequencies unless told otherwise
if (sides == 'default' and np.iscomplexobj(x)) or sides == 'twosided':
num_freqs = pad_to
scaling_factor = 1.
elif sides in ('default', 'onesided'):
num_freqs = pad_to // 2 + 1
scaling_factor = 2.
else:
raise ValueError("sides must be one of: 'default', 'onesided', or "
"'twosided'")
# Matlab divides by the sampling frequency so that density function
# has units of dB/Hz and can be integrated by the plotted frequency
# values. Perform the same scaling here.
if scale_by_freq:
scaling_factor /= Fs
window_vals = np.hanning(NFFT)
step = int(NFFT) - int(noverlap)
ind = np.arange(0, len(x) - NFFT + 1, step, dtype=np.int32)
n = len(ind)
p_xy = np.zeros((num_freqs, n), np.complex_)
# do the ffts of the slices
for i in range(n):
this_x = x[ind[i]:ind[i] + NFFT]
this_x = window_vals * this_x
fx = np.fft.fft(this_x, n=pad_to)
if same_data:
fy = fx
else:
th_is_y = y[ind[i]:ind[i] + NFFT]
th_is_y = window_vals * th_is_y
fy = np.fft.fft(th_is_y, n=pad_to)
p_xy[:, i] = np.conjugate(fx[:num_freqs]) * fy[:num_freqs]
# Scale the spectrum by the norm of the window to compensate for
# windowing loss; see Bendat & Piersol Sec 11.5.2. Also include
# scaling factors for one-sided densities and dividing by the sampling
# frequency, if desired.
p_xy *= scaling_factor / (np.abs(window_vals) ** 2).sum()
t = 1. / Fs * (ind + NFFT / 2.)
freqs = float(Fs) / pad_to * np.arange(num_freqs)
if (np.iscomplexobj(x) and sides == 'default') or sides == 'twosided':
# center the frequency range at zero
freqs = np.concatenate((freqs[num_freqs // 2:] - Fs,
freqs[:num_freqs // 2]))
p_xy = np.concatenate((p_xy[num_freqs // 2:, :],
p_xy[:num_freqs // 2, :]), 0)
return p_xy, freqs, t