r"""

# Use numpy to find where thresholds are exceeded

mc_good = np.where(mdccm > mcthresh)[0]

az_good = np.where((az >= az_volc - az_diff) & (az <= az_volc + az_diff))[0]

vel_good = np.where((vel >= .25) & (vel <= .45))[0]

igood = reduce(np.intersect1d, (mc_good, az_good, vel_good))

# Find find number of consecutive values exceeded.

ranges = []

nconsec = []

for k, g in groupby(enumerate(igood), lambda x: x[0]-x[1]):

group = list(map(itemgetter(1), g))

ranges.append((group[0], group[-1]))

nconsec.append(group[-1]-group[0]+1)

if len(nconsec) > 0:

consecmax = max(nconsec)

else:

consecmax = 0

print('%d above trheshold, %d consecutive\n' % (len(igood), consecmax))

M=None, Moffset=None):

r"""

# (c) 2017 Curt A. L. Szuberla

# University of Alaska Fairbanks, all rights reserved

#

# size things up

m, nTraces = data.shape

# -- default parsing & error checking -----------------------------------

# default weight is unity weight for all traces

if wgt is None:

wgt = [1 for i in range(nTraces)]

else:

if len(wgt) != nTraces:

# catch dimension mismatch between tau & wgt

raise IndexError('len(wgt) != ' + str(nTraces))

wgt = np.array(wgt) # require array form here for later operations

# default refTrace is first non-zero wgt

if refTrace is None:

refTrace = np.min(np.where(wgt != 0)) # requires array wgt

# default Moffset is zero for all traces

if Moffset is None:

Moffset = [0 for i in range(nTraces)]

else:

if len(Moffset) != nTraces:

# catch dimension mismatch between tau & Moffset

raise IndexError('len(Moffset) != ' + str(nTraces))

# -- end default parsing & error checking -------------------------------

# planar (far-field) or spherical (near-field) arrival?

if r is None:

tau = tauCalcPW(vel, azPhi, rij)

else:

# need to unpack & repack azPhi with care

if np.isscalar(azPhi):

tau = tauCalcSW(vel, [r, azPhi], rij)

else:

tau = tauCalcSW(vel, [r, azPhi[0], azPhi[1]], rij)

# calculate shifts as samples

beam_delays = phaseAlignIdx(tau, Hz, wgt, refTrace)

# apply shifts, resulting in a zero-padded array

beamMatrix = phaseAlignData(data, beam_delays, wgt, refTrace, M, Moffset)

# linear algrebra to perform sum & then normalize by weights

r"""

# (c) 2017 Curt A. L. Szuberla

# University of Alaska Fairbanks, all rights reserved

#

# -- this is low level code w/out error checks or defaults, designed

# -- to be called by wrappers that make use of the indices provided

# size up data

m, nTraces = data.shape

# if plotting, embed in array of np.nan

if plotFlag:

nanOrOne = np.nan

else:

nanOrOne = 1

# correct for negative Moffset elements

# subtract this to ensure corrected delays is positive,

# semi-definite and has (at least) one zero element

maxNegMoffset = min(np.array(Moffset)[np.array(Moffset) <= 0])

# apply Moffset & correction for negative elements of Moffset

delays = delays + Moffset - maxNegMoffset

# start with everything in the window as a default (trim||pad later)

data_align = np.zeros((max(delays) + m, nTraces)) * nanOrOne

# embed shifted traces in array

for k in range(nTraces):

if wgt[k]:

data_align[delays[k]:delays[k]+m, k] = data[:, k] * wgt[k]

# truncate|| pad data_align if M >< m, centered on refTrace

mp = data_align.shape[0] # new value for m

if M is not None and M is not mp:

alignBounds = [delays[refTrace] + m//2 - M//2,

delays[refTrace] + m//2 + M//2]

# trap round-off errors and force (M, nTraces) data_align

if alignBounds[1]-alignBounds[0] != M:

alignBounds[1] += 1

if not (alignBounds[1]-alignBounds[0])%2:

alignBounds[0] -= 1

# -- LHS (graphically, but actually topside in array-land!)

if alignBounds[0] < 0:

# pad LHS of traces w zeros or np.nans

data_align = np.vstack((np.zeros((-alignBounds[0], nTraces)) * nanOrOne,

data_align))

elif alignBounds[0] > 0:

data_align = data_align[alignBounds[0]:]

# -- RHS (graphically, but actually bottom in array-land!)

if alignBounds[1] > mp:

# pad RHS of traces w zeros or np.nans

data_align = np.vstack((data_align, np.zeros((alignBounds[1] - mp,

nTraces)) * nanOrOne))

elif alignBounds[1] < mp:

data_align = data_align[:M]

r"""

# -- this is low level code w/out error checks or defaults, designed

# -- to be called by wrappers that make use of the indices provided

# solve for number of traces from pairings in tau

nTraces = int(1 + np.sqrt(1 + 8 * len(tau))) // 2

# calculate delays (samples) relative to refTrace for each trace

# -- first pass grabs delays starting with refTrace as i in ij

delayIdx = (nTraces*refTrace - refTrace*(refTrace+1)//2,

nTraces*(refTrace+1) - (refTrace+1)*(refTrace+2)//2)

delays = np.hstack((0, (tau[range(delayIdx[0], delayIdx[1])] * Hz))).astype(int)

# the std. rij list comprehension for unique inter-trace pairs

tau_ij = [(i, j) for i in range(nTraces) for j in range(i+1, nTraces)]

# -- second pass grabs delays with refTrace as j in ij

preRefTau_idx = [k for k in range(len(tau)) if tau_ij[k][1] == refTrace]

delays = np.hstack((-tau[preRefTau_idx] * Hz, delays)).astype(int)

# re-shift delays such that the first sample of the trace requiring the

# largest shift left relative to the refTrace (hence, largest positive,

# semi-definite element of delays) has an index of zero; i.e., all traces

# have a valid starting index for embedding into an array (no padding)

r"""

dim, nTraces = rij.shape

if dim == 2:

rij = np.vstack((rij, np.zeros((1, nTraces))))

idx = [(i, j) for i in range(rij.shape[1]-1)

for j in range(i+1, rij.shape[1])]

X = rij[:, [i[0] for i in idx]] - rij[:, [j[1] for j in idx]]

if np.isscalar(azPhi):

phi = 0

az = azPhi

else:

phi = azPhi[1] / 180 * np.pi

az = azPhi[0]

az = np.pi * (-az / 180 + 0.5)

s = np.array([np.cos(az), np.sin(az), np.sin(phi)])

s[:-1] *= np.cos(phi)

r"""

dim, nTraces = rij.shape

if dim == 2:

rij = np.vstack((rij, np.zeros((1, nTraces))))

if len(rAzPhi) == 3:

phi = rAzPhi[2] / 180 * np.pi

else:

phi = 0

idx = [(i, j) for i in range(rij.shape[1]-1)

for j in range(i+1, rij.shape[1])]

# aw, this is so convolutedly elegant that it must be saved in a

# comment for posterity!, but the line below it is "simpler"

# az = -( (rAzPhi[1]/180*pi - 2*pi)%(2*pi) - pi/2 )%(2*pi)

az = np.pi * (-rAzPhi[1] / 180 + 0.5)

source = rAzPhi[0] * np.array([np.cos(az), np.sin(az), np.sin(phi)])

source[:-1] *= np.cos(phi)

tau2sensor = np.linalg.norm(rij - np.tile(source, nTraces).reshape(nTraces, 3).T,

2, axis=0)/vel

r"""

dim, nTraces = len(rij), len(rij[0])

if dim == 2:

rij = np.vstack((rij, [0] * nTraces))

else:

rij = np.vstack((rij,))

if len(xy) == 2:

xy0 = 0

else:

xy0 = []

source = np.hstack((xy, xy0))

idx = [(i, j) for i in range(rij.shape[1]-1)

for j in range(i+1, rij.shape[1])]

tau2sensor = np.linalg.norm(rij - np.tile(source, nTraces).reshape(nTraces, 3).T,

2, axis=0)/vel

r"""

# catch scalar input & form tuple

if type(N) is int:

N = (N,)

# ensure DC component of output will be zero (unless beta == 0, see below)

if beta < 0:

c0 = 0

else:

c0 = np.inf

# catch the case of a 1D array in python, so dimensions act like a matrix

if len(N) == 1:

M = (N[0], 1) # use M[1] any time # of columns is called for

else:

M = N

# phase array with size (# of unique complex Fourier components,

# columns of original data)

n = int(np.floor((N[0] - 1) / 2)) # works for odd/even cases

cPhase = np.random.random_sample((n, M[1])) * 2 * np.pi

# Nyquist placeholders

if N[0]%2:

# odd case: Nyquist is 1/2 freq step between highest components

# so it is empty

cFiller = np.empty((0,))

pFiller = np.empty((0, M[1]))

else:

# even case: we have a Nyquist component

cFiller = N[0]/2

pFiller = np.zeros((1, M[1]))

# noise amplitudes are just indices (unit-offset!!) to be normalized

# later, phases are arranged as Fourier conjugates

r = np.hstack((c0, np.arange(1, n + 1), cFiller, np.arange(n, 0, -1)))

phasor = np.exp(np.vstack((np.zeros((1, M[1])), 1j * cPhase, pFiller,

-1j * np.flipud(cPhase))))

# this is like my cols.m function in MATLAB

r = np.tile(r, M[1]).reshape(M[1], N[0]).T ** (-beta / 2)

# catch beta = 0 case here to ensure zero DC component

if not beta:

r[0] = 0

# inverse transform to get time series as columns, ensuring real output

X = r*phasor

r = np.real(np.fft.ifft(X, axis=0)*X.shape[0])

# renormalize r such that mean = 0 & std = 1 (MATLAB dof default used)

# and return it in its original shape (i.e., a 1D vector, if req'd)

r"""

# private functions up front

def Ssmooth(S, w, n, window):

# smooth special format of spectral matries as vectors

for k in range(n):

# f@#$ing MATLAB treats odd/even differently with mode='full'

# but the behavior below now matches conv2 exactly

S = convolve2d(S, window(w).reshape(-1, 1),

mode='full')[w//2:-w//2+1, :]

return S

def triang(N):

# for historical reasons, the default window shape

return np.bartlett(N + 2)[1:-1]

# size up input data

N, d = x.shape

# Fourier transform of data matrix by time series columns, retain only

# the diagonal & above (unique spectral components)

Nx = x.shape

X = np.fft.fft(x, axis=0)/Nx[0]

X = X[:N//2+1, :]

# form spectral matrix stack in reduced vector form (**significant**

# speed improvement due to memory problem swith 3D tensor format -- what

# was too slow in 1995 is still too slow in 2017!)

# -- pre-allocate stack & fill it in

Sidx = [(i, j) for i in range(d) for j in range(i, d)]

S = np.empty((X.shape[0], d * (d + 1) // 2), dtype=complex)

for i in range(X.shape[0]):

# at each freq w, S_w is outer product of raw Fourier transforms

# of each column at that freq; select unique components to store

S_w = np.outer(X[i, :], X[i, :].conj())

S[i, :] = S_w[[i[0] for i in Sidx], [j[1] for j in Sidx]]

# smooth each column of S (i,e., in freq. domain)

if not window:

# use default window

window = triang

S = Ssmooth(S, w, n, window)

# trace calculations (notation consistent w traceCalc.m in MATLAB), but

# since results are positive, semi-definite, real -- return as such

# -- diagonal elements

didx = [i for i in range(len(Sidx)) if Sidx[i][0] == Sidx[i][1]]

# -- traceS**2 of each flapjack (really a vector here)

trS = sum(S[:, didx].real.T)**2

# -- trace of each flapjack (ditto, vector), here we recognize that

# trace([email protected]) is just sum square magnitudes of all the

# non-redundant components of S, doubling squares of the non-diagonal

# elements

S = (S*(S.conj())*2).real

S[:, didx] /= 2

trS2 = sum(S.T)

# estimate samson-esque polarization estimate (if d==2, same as fowler)

P = (d*trS2 - trS)/((d-1)*trS)

# a litle trick here to handle odd/even number of samples and zero-out

# both the DC & Nyquist (they're both complex-contaminated due to Ssmooth)

P[0] = 0

if N%2:

# odd case: Nyquist is 1/2 freq step between highest components

fudgeIdx = 0

else:

# even case: we have a Nyquist component

fudgeIdx = 1

P[-1] = 0

# apply P as contrast agent to frequency series

X *= np.tile(P ** p, d).reshape(X.shape[::-1]).T

# inverse transform X and ensure real output

XX = np.vstack((X[list(range(N // 2 + 1))],

X[list(range(N//2-fudgeIdx, 0, -1))].conj()))

x_psf = np.real(np.fft.ifft(XX, axis=0)*XX.shape[0])