array_processing.tools.generic module

array_processing.tools.generic.array_thresh(mcthresh, az_volc, az_diff, mdccm, az, vel)[source]

Find array processing values above multiple set thresholds for MCCM, back-azimuth, and trace velocity. Uses default 0.25–0.45 km/s for trace velocity thresholds. Also finds consecutive segments that meet thresholds, but these values are not currently returned.

Parameters:
  • mcthresh (float) – MCCM or MdCCM threshold (0–1)
  • az_volc (float) – Back-azimuth to target volcano or source (0–359)
  • az_diff (float) – Tolerance for back-azimuth from az_volc
  • mdccm – MdCCM or MCCM values from array processing (0–1)
  • az – Back-azimuth values (0–359)
  • vel – Trace-velocity values (km/s)
Returns:

Indices to time segments that meet set thresholds
array_processing.tools.generic.beamForm(data, rij, Hz, azPhi, vel=0.34, r=None, wgt=None, refTrace=None, M=None, Moffset=None)[source]

Form a “best beam” from the traces of an array.

Parameters:
  • data – (m, n) array; time series with m samples from n traces as columns
  • rij – (d, n) array; n sensor coordinates as [easting, northing, {elevation}] column vectors in d dimensions
  • Hz (int or float) – Sample rate
  • azPhi – Back azimuth (float) from co-array coordinate origin (° CW from N); back azimuth and elevation angle (list) from co-array coordinate origin (° CW from N, ° from N-E plane)
  • vel (float) – Estimated signal velocity across array
  • r (float) – Range to source from co-array origin. Default is None (use plane wave arrival model), If not None, use spherical wave arrival model
  • wgt – Vector of relative weights of length n (0 == exclude trace). Default is None (use all traces with equal relative weights [1 for i in range(nTraces)])
  • refTrace (int) – Reference sensor for TDOA information. Default is None (use first non-zero-weighted trace)
  • M (int) – Length of best beam result in samples. Default is None (use m samples, same as input data)
  • Moffset – Individual trace offsets from arrival model shifts. Default is None (use [0 for i in range(nTraces)])
Returns:

(M, ) array of summed and weighted shifted traces to form a best beam
Raises:

IndexError – If the input argument dimensions are not consistent

Notes

This beamformer handles planar- or spherical-model arrivals from arbitrarily elevated sources incident on 2- or 3-D arrays. Weights are relative and normalized upon beam output. The default value for vel assumes that rij is in units of km (e.g., the speed is in km/s).

array_processing.tools.generic.phaseAlignData(data, delays, wgt, refTrace, M, Moffset, plotFlag=False)[source]

Embeds n phase aligned traces in a data matrix.

Parameters:
  • data – (m, n) array; time series with m samples from n traces as columns
  • delays – (n, ) array; vector of shifts as indices for embedding traces in an array, such that trace i will begin at index out[i]
  • wgt – Vector of relative weights of length n (0 == exclude trace by setting to padding value, see plotFlag)
  • refTrace (int) – Reference sensor for TDOA information
  • M (int) – Length of best beam result in samples (use m to let beam be same length as input traces)
  • Moffset – Individual trace offsets from arrival model shifts (use [0 for i in range(nTraces)] to skip this effect)
  • plotFlag (bool) – Flag to indicate output array will be used for plotting purposes. Default is False (pads shifts with zeros; pads with numpy.nan if True)
Returns:

(M, n) array of shifted traces as columns

Notes

The output of phaseAlignIdx() is used to calculate the input delays.

array_processing.tools.generic.phaseAlignIdx(tau, Hz, wgt, refTrace)[source]

Calculate shifts required to phase align n traces in a data matrix.

Parameters:
  • tau – (n(n-1)//2, ) array; time delays of relative signal arrivals (TDOA) for all unique sensor pairings
  • Hz (int or float) – Sample rate
  • wgt – Vector of relative weights of length n (0 = exclude trace)
  • refTrace (int) – Reference sensor for TDOA information
Returns:

(n, ) array; vector of shifts as indices for embedding traces in an array, such that trace i will begin at index out[i]

Notes

The output of this function is compatible with the inputs of phaseAlignData().

array_processing.tools.generic.psf(x, p=2.0, w=3, n=3.0, window=None)[source]

Pure-state filter a data matrix. This function uses a generalized coherence estimator to enhance coherent channels in an ensemble of time series.

Parameters:
  • x – (m, d) array of real-valued time series data as columns
  • p (float) – Level of contrast in filter
  • w (int) – Width of smoothing window in frequency domain
  • n (float) – Number of times to smooth in frequency domain
  • window – Type of smoothing window in frequency domain. Default is None, which results in a triangular window
Returns:

Tuple containing:

  • x_psf – (m, d) array; real-valued, pure state-filtered version of x
  • P – (m//2+1, ) array; degree of polarization (generalized coherence estimate) in frequency components of x from DC to the Nyquist
Return type:

tuple

Notes

See any of Samson & Olson’s early 1980s papers, or Szuberla’s 1997 PhD thesis, for a full description of the underlying theory. The code implementation’s defaults reflect historical values for the smoothing window — a more realistic w would be of order \(\sqrt{m}\) combined with a smoother window, such as numpy.hanning(). Letting n=3 is a reasonable choice for all window types to ensure confidence in the spectral estimates used to construct the filter.

For \(m\) samples of \(d\) channels of data, a (d, d) spectral matrix \(\mathbf{S}[f]\) can be formed at each of the m//2+1 real frequency components from DC to the Nyquist. The generalized coherence among all of the \(d\) channels at each frequency is estimated by

\[P[f] = \frac{d \left(\text{Tr}\,\mathbf{S}^2[f]\right) - \left(\text{Tr}\,\mathbf{S}[f]\right)^2} {\left(d-1\right)\left(\text{Tr}\,\mathbf{S}[f]\right)^2},\]

where \(\text{Tr}\,\mathbf{S}[f]\) is the trace of the spectral matrix at frequency \(f\). The filter is constructed by applying the following multiplication in the frequency domain

\[\hat{\mathbf{X}}[f] = P[f]^p\mathbf{X}[f],\]

where \(\mathbf{X}[f]\) is the Fourier transform component of the all channels at frequency \(f\) and \(p\) is the level of contrast. The inverse Fourier transform of the matrix \(\hat{\mathbf{X}}\) gives the filtered time series.

The estimator \(\mathbf{P}[f] = 1\), identically, without smoothing in the spectral domain (a consequence of the variance in the raw Fourier components), but it is bound by \(\mathbf{P}[f]\in[0,1]\) even withn smoothing, hence its utility as a multiplicative filter in the frequency domain. Similarly, this bound allows the contrast between channels to be enhanced based on their generalized coherence if \(p>1\).

Data channels should be pre-processed to have unit-variance, since unlike the traditional two-channel magnitude squared coherence estimators, the generalized coherence estimate can be biased by relative amplitude variations among the channels. To mitigate the effects of smoothing complex values into the DC and Nyquist components, they are set to zero before computing the inverse transform of \(\hat{\mathbf{X}}\).

array_processing.tools.generic.randc(N, beta=0.0)[source]

Colored noise generator. This function generates pseudo-random colored noise (power spectrum proportional to a power of frequency) via fast Fourier inversion of the appropriate amplitudes and complex phases.

Parameters:
  • N (int or tuple) – Shape of output array
  • beta (float) – Spectrum of output will be proportional to f**(-beta)
Returns:

Colored noise sequences as columns with shape N, each normalized to zero-mean and unit-variance. Result is always real valued

Notes

Spectrum of output will be \(\sim1/f^\beta\). White noise is the default (\(\beta = 0\)); others are pink (\(\beta = 1\)) or brown/surf (\(\beta = 2\)). Since the output is zero-mean, the DC spectral component(s) will be identically zero.

array_processing.tools.generic.tauCalcPW(vel, azPhi, rij)[source]

Calculates theoretical tau vector for a plane wave moving across an array of n elements.

Parameters:
  • vel (float) – Signal velocity across array
  • azPhi – Back azimuth (float) from co-array coordinate origin (° CW from N); back azimuth and elevation angle (array) from co-array coordinate origin (° CW from N, ° from N-E plane)
  • rij – (d, n) array; n element coordinates as [easting, northing, {elevation}] column vectors in d dimensions
Returns:

(n(n-1)//2, ) array; time delays of relative signal arrivals (TDOA) for all unique sensor pairings
array_processing.tools.generic.tauCalcSW(vel, rAzPhi, rij)[source]

Calculates theoretical tau vector for a spherical wave moving across an array of n elements.

Parameters:
  • vel (float) – Signal velocity across array
  • rAzPhi – Range to source and back azimuth from co-array coordinate origin (° CW from N); range to source, back azimuth and elevation angle from co-array coordinate origin (° CW from N, ° from N-E plane)
  • rij – (d, n) array; n element coordinates as [easting, northing, {elevation}] column vectors in d dimensions
Returns:

(n(n-1)//2, ) array; time delays of relative signal arrivals (TDOA) for all unique sensor pairings
array_processing.tools.generic.tauCalcSWxy(vel, xy, rij)[source]

Calculates theoretical tau vector for a spherical wave moving across an array of n elements.

Parameters:
  • vel (float) – Signal velocity across array
  • xy – (d, ) array; source location as 2-D [easting, northing] or 3-D [easting, northing, elevation] coordinates
  • rij – (d, n) array; n element coordinates as [easting, northing, {elevation}] column vectors in d dimensions
Returns:

(n(n-1)//2, ) array; time delays of relative signal arrivals (TDOA) for all unique sensor pairings