coodddaaaa.stretching

Copyright (c) 2023 maximilien.lehujeur

class coodddaaaa.stretching.Stretcher(t0: float, dt: float, nt: int, eps: ndarray, norm: bool = False, interp_kind: Literal['linear', 'cubic', 'fourier'] = 'cubic')

Bases: object

An object to compute stretched signals and to perform stretching correlation as defined in Weaver et al., 2011.

\[X(\varepsilon) = \frac {\int{y^{ref}(t \times (1 + \varepsilon)) \cdot y(t) dt}} {\sqrt{ \int{y^{ref}(t \times (1 + \varepsilon))^2 dt} \int{y(t)^2 dt} }}\]

The object pre-computes the interpolation operator. The user can pre-compute and store the stretched basis of the reference signal. This basis can then be provided for stretching correlation with a new signal. This object can also compute the stretching between all pairs of signals in a b-scan.

Parameters:

• t0 – time of first sample

• dt – sampling interval

• nt – number of samples

• eps – epsilon array

• norm – use it to compute normalized correlation. warning : for stretching only, use norm = False

• interp_kind – which interpolator to use for stretching, among ‘linear’, “cubic”, ‘fourier’

stretch(x: ndarray) → ndarray

Compute the stretched basis functions from a signal x

Parameters:

x – the input signal (reference), np.ndarray, 1d, shape (nt, )

Return x_stretched:

the stretched version of x for all values in self.eps, np.ndarray 2d, shape (neps, nt)

corr(x: ndarray, x_stretched: ndarray) → ndarray

Stretching correlation of x with a basis of stretched versions of the reference signal

Parameters:

• x – signal(s) to be correlated to the reference, np.ndarray, either one single signal, 1d, shape (nt, ) or a bscan, 2d, shape (ntraces, nt)

• x_stretched – stretched reference from self.stretch, np.ndarray, 2d, shape (neps, nt, )

Return c:

correlation function np.ndarray, either 1d, shape (neps, ) if x is 1d or 2d, shape (neps, ntraces) if x is 2d

corrmax(c: ndarray)

Find the maximum of the correlation function with subsample precision

Parameters:

c – correlation function(s) from self.corr 1d for a single signal, shape (neps, ) 2d for a bscan, shape (neps, ntraces)

Return emax:

best epsilon value, dimensionless, it corresponds to dt/t float if c is 1d 1d array, shape (ntraces, ) if c is 2d

Return cmax:

max correlation, dimensionless, normalized if norm was True in __init__ float if c is 1d 1d array, shape (ntraces, ) if c is 2d

corr_all_with_all(data: ~numpy.ndarray) -> (<class 'numpy.ndarray'>, <class 'numpy.ndarray'>)

Correlate all possible pairs of signals in a bscan

Parameters:

data – the bscan, one trace per row, same sampling (=self.t), 2d, shape (ntraces, nt)

Return c_triu:

the max correlation coefficients for all pairs (upper triangle only)

Return e_triu:

the best stretching coefficients for all pairs (upper triangle only)

use self.triu2dence to get the full matrices

c = Stretcher.triu2dense(c_triu, symetric=True, diag=1.0) e = Stretcher.triu2dense(e_triu, symetric=False, diag=0.0)

static triu2dense(x_triu: ndarray, symetric: bool, diag: float) → ndarray

Convert upper triangle matrix to square matrix

Parameters:

• x_triu – a flat upper triangle without diagonal, 1d, np.ndarray, shape (ntraces * (ntraces - 1) / 2, )

• symetric – to impose symetry (True) or anti-symetry (False)

• diag – the value to put on the diagonal

Return x:

a square matrix with x_triu on its upper triangle, shape (ntraces, ntraces) diag on its diagonal +-x_triu on its lower triangle

static stretching_uncertainty(cmax: float | ndarray, fmin: float, fmax: float, tmin: float, tmax: float) → float | ndarray

Stretching uncertainty after Weaver et al 2011.

Parameters:

• cmax – max correlation coefficient from self.corrmax, either a float or a np.ndarray

• fmin – lower freq Hz, float

• fmax – upper freq Hz, float

• tmin – start coda time in s, float

• tmax – end coda time in s, float

Return rmse:

uncertainty on epsilon, same type as cmax

class coodddaaaa.stretching.InverseStretcher(t0: float, nt: int, dt: float, eps_history: ndarray, interp_kind: Literal['linear'] = 'linear')

Bases: object

An object to cancel the effect of the stretching on each trace of a bscan This can be used to align the traces with the reference, and then to refine the reference.

For example:

you have a bscan of 256 traces with n samples each, bscan is a 2d array shapped (256, n) you have an estimate of the stretching history, i.e. 256 epsilon values in an 1D array this object returns the bscan corrected from the estimated stretching values

positive epsilon values (i.e. positive dv/v) mean that the trace was compressed relative to its ref, so this operator stretch it negative epsilon values will tend to compress the waveform

Parameters:

• t0 – time of first sample

• nt – number of samples

• dt – sampling interval

• eps_history – epsilon array, one item per trace in the bscan

• interp_kind – which interpolator to use for inverse stretching, among ‘linear’,

coodddaaaa.interp1d

Copyright (c) 2023 maximilien.lehujeur

Linear and cubic interpolation in 1d using fixed grids and sparse operators for cases where one need to interpolate functions on the same grids many times

note I do not use the scipy interpolator because I need an interpolator that can be created from the grids only and called later on with the function to interpolate

2023.04.07 : P. Mora : Speed up the construction of the sparse matrixes => x20 to x50

class coodddaaaa.interp1d.LinearInterpolator1d(x0: float, nx: int, dx: float, xi: ndarray, format: str = 'csc')

Bases: object

Linear interpolation operator

Parameters:

• x0 – x coordinate of the first sample

• nx – number of samples

• dx – sampling interval

• xi – the points where we need the interpolated values => f(xi) is computed by self.__call__

• format – format to use for the linear operator

find_interp_points_in_grid(x0: float, nx: int, dx: float, xi: ndarray)

class coodddaaaa.interp1d.SecondDerivativeOperatorTypeII(nx: int, dx: float, format: str = 'csc')

Bases: object

Second derivative operator order 3 in the internal domain, Implement type II boundary condition after https://en.wikiversity.org/wiki/Cubic_Spline_Interpolation For a regular grid only x is the grid at which the function will be defined (nodes) xi are the points where the function will be interpolated

Parameters:

• nx – number of nodes

• dx – sampling interval between nodes

• format – format of the sparse operator

class coodddaaaa.interp1d.CubicInterpolator1d(x0: float, nx: int, dx: float, xi: ndarray, format: str = 'csc')

Bases: LinearInterpolator1d

Lagrange Cubic interpolation with boundary type II from https://en.wikiversity.org/wiki/Cubic_Spline_Interpolation Works only on a regular grid for now x is the grid at which the function will be defined (nodes) xi are the points where the function will be interpolated :param x0: x of first sample :param nx: number of samples :param dx: sampling interval :param xi: array of points where to interpolate the function :param format: format of the sparse operator

class coodddaaaa.interp1d.RFFTInterpolator1d(x0: float, nx: int, dx: float, xi: ndarray)

Bases: object

Fourier Interpolator based on rfft

Parameters:

• x0 – x of first sample

• nx – number of samples

• dx – sampling interval

• xi – array of points where to interpolate the function

• format – format of the sparse operator

coodddaaaa.utils

Copyright (c) 2023 maximilien.lehujeur

class coodddaaaa.utils.Timer(message: str)

Bases: object

Counts the execution time under the “with” statement

Parameters:

message

coodddaaaa.utils.polyspace(xmin: float, xmax: float, nx: int, pwr: float)

A power-law to refine resolution of a stretching grid search near zero

Parameters:

• xmin – min value

• xmax – max value

• nx – number of points

• pwr – power coefficient, increase pwr to refine the resolution near 0

class coodddaaaa.utils.TukeyWindow(t0: float, t1: float, t2: float, t3: float)

Bases: object

A parameterizable 4 points Tukey function :param t0 … t3: times of the corners of the Tukey window

coodddaaaa.butter

Modified after sigy 1.5.3, M.L. 21/04/2023

Time / Fourier domain butterworth filter warning : the fourier domain filter has a slightly different response that the time domain one

comparing both reveals that the time domain filter may include a water-level that do not not exist with fourier domain, this results in slight differences near the signal edges taper the waveform properly fixes the difference

TODO : use ba_analog ? TODO : use scipy.fft => parallel

class coodddaaaa.butter.ButterworthFilter(freqmin: float | None, freqmax: float | None, sampling_rate: float | None, order: float = 4.0)

Bases: object

Parameters:

• freqmin – lower frequency in Hz, or None for highpass filtering

• freqmax – upper frequency in Hz, or None for lowpass filtering

• sampling_rate – in Hz

• order – of the filter

timecall(data, zerophase=False, axis=-1)

response(npts, zerophase=False, input_domain='fft', qc=False)

Almost equivalent to timecall the response looks better than timecall, no water level applied

show(fig, freqs=None, zerophase=False, **kwargs)

class coodddaaaa.butter.BandpassFilter(freqmin, freqmax, sampling_rate, order=4)

Bases: ButterworthFilter

Shortcut for ButterworthFilter for band-pass filtering

class coodddaaaa.butter.LowpassFilter(freqmax, sampling_rate, order=4)

Bases: ButterworthFilter

Shortcut for ButterworthFilter for low-pass filtering

class coodddaaaa.butter.HighpassFilter(freqmin, sampling_rate, order=4)

Bases: ButterworthFilter

Shortcut for ButterworthFilter for high-pass filtering

coodddaaaa.fftoversamp

Modified after sigy 1.5.3, M.L. 20/08/2023

Fourier domain oversampling, for the sake of simplicity, this program can increase the number of samples only by 2**n, where n is an integer.

The oversampling is performed in the Fourier domain, - for optimal use, it is preferred to use the fft_oversamp

if the signal is already in Fourier domain

• make sure the signal is properly detrended and tapered at its edges prior to FFT otherwise, you might observe wiggles at the edges of the signal after oversampling. (remember that fft assume a periodization of the signal in time)

coodddaaaa.fftoversamp.oversamp(t0: float, dt: float, data: ~numpy.ndarray, npow2: int, axis: int = -1, demean: bool = False) -> (<class 'numpy.ndarray'>, <class 'numpy.ndarray'>)

Time domain version of fft_oversamp

Parameters:

• t0 – start time, sec

• dt – sampling interval, sec

• data – time domain data array (1d or more)

• npow2 – oversamp by 2 ** npow2

• axis – axis along which to oversample the signal

Return to:

the new time vector

Return datao:

the oversample data

coodddaaaa.fftoversamp.fft_oversamp(fft_data: ndarray, npow2: int = 1, axis: int = -1) → ndarray

Oversamp a signal by padding it with zeros in the FFT domain the number of sample is multiplied by 2 ** npow2 (default 2**1)

Parameters:

• fft_data – output of fft

• npow2 – oversampling rate expressed as a power of 2

• axis – the axis along which to perform oversampling

Returns:

oversample data in fft domain

coodddaaaa.hypermax

Copyright (c) 2023 maximilien.lehujeur

Finds the maximum of a random array with subsample precision by looking for the zero crossing of the first order finite difference derivative

coodddaaaa.hypermax.hypermax(time_array: ndarray, function_array: ndarray, axqc: Axes | None = None, assume_t_growing: bool = False)

Parameters:

• time_array – time array

• function_array – function array

• axqc – matplotlib ax or None for visual qc

• assume_t_growing