xcor module

This optional module allows to compute cross-correlation from a 2D ndarray (or a A1Section) containing a DAS section. The core function is written in fortran and parallelized using openMP.

Instead of processing cross-correlation using individual records, the xcor module uses directly the 2D float array containing DAS data. A pre-processing step to select the couple of traces to be cross-correlated is achieved by calling the register_couple() function. A number of pre-processing are built-in and available (one-bit, whitenning, bandpass filtering, …) by calling the register_par() function.

Requirements

This module requires:

  • C and fortran compilers (e.g gfortran) and make (gmake is ok)

  • numpy

  • one of the MKL or OpenBlas+FFTW3 libraries.

For linux, these are standard and free libraries that can be easily download using apt-get. Try; apt-get install libfftw3 apt-get install libfftw3-dev apt-get install libopenblas-dev

For mac with arm cpu (M1, M2, …) using OpenBlas and FFTW is mandatory, they can be get from brew. Getting a fortran compiler usually requires to install Apple xcode.

Download and installation

  1. Install (tar.gz archive) the xcor module in the same directory as a1das from the gitlab repository (e.g. in2p3 gitlab), you should see a directory named A1-XCORPY.

  2. Check that the required libraries are installed and find their locations (/opt/…, /usr/local/lib, /usr/lib/x86_64-linux-gnu, /opt/homebrew, ….). You should look for libopenblas.a, libfftw3.a or libmkl

  3. change to A1-XCORPY directory and under a python environment, run the shell configure_var.sh, this will create the Makefile.python file

  4. select one of the proposed Makefile.inc.???, edit and change according to your needs, then copy it as Makefile.inc

  5. run the compilation step by typing “make python”

This should have installed the binary module under the a1das directory. If you have already installed the a1das package you can either

a) uninstall and reinstall it: python -m pip uninstall a1das pip install .

b) or copy the _a1xcorPy.cpython.?? file in the directory where a1das has been installed (e.g.

Usage

xcor.compute_xcorr(arrayIn, lag=None, stack=False, verbose=0, base=0, flip=False)

Description

xcorr, lags, ijx, ier = compute_xcorr(arrayIn, lag=None, stack=False, verbose=0)

Compute cross-correlations in a 2D array of traces [Ntime x Nspace] The array is given either as a 2D ndarray or as an A1Section class instance.

Cross-correlations are computed for the couple of traces {i,j} declared previously by a1das.xcor.register_couple() and using the parameter declared previously by a1das.xcor.register_par()

!!!!! WARNING !!!! currently process only float64 data array

TODO: write a float32 version

Input

arrayIn:

2D (float64) ndarray [nSpace x nTime], or <A1Section> class instance (with float64 data array)

lag:

(float) time lag, in sample for ndarray input, in sec for <A1Section>

stack:

(bool) stack the correlations for consecutive windows of 2*lag+1 length if 2*lag+1 < length signal (default False)

base:

(int) <0> if first index of trace is 0, <1> if first index is 1 (default 0)

flip:

(bool) reverse time lags with respect to usual python correlation (default false)

verbose:

verbosity (default 0)

function return

xcorr:

2D ndarray of cross-correlations [nxcor x nTime] (float64)

lags:

vector of lag times, in sample for ndarray input, in sec for <A1Section> input

ijx:

2D ndarray of size[nxcor, 2] giving trace index for each correlation

ier:

return code, 0 if ok

Usage example

>>> import a1das
>>> f = a1das.open(filename, format='febus')
>>> a1 = f.read()
>>> a1das.xcor.register_couple([[1,2],[1,3],[1,4]]) # compute xcorr between traces 1&2, 1&3, 1&4
>>> a1das.xcor.register_par('onebit') # set onebit pre-processing
>>> lag = 5.           # set lag time to 5sec
>>> # compute xcorr between -5sec and 5sec, compute FFT over consecutive windows of 10 sec and stack them
>>> xcorr, lags, ijx, ier = compute_xcorr(a1, lag=lag, stack = True)
>>> # same but for ndarray with no stack,
>>> xcorr, lags, ijx, ier = compute_xcorr(a1.data, lag=5./a1['dt'], stack = False)
xcor.register_couple(couple, ntrace, base=0)

Description

Register the list of trace couples {i,j} for which we compute the cross-correlation

if base = 0 (default = python, C convention)

0 <= i < ntrace 0 <= j < ntrace

if base = 1 (Fortran, matlab convention)

1 <= i <= ntrace 1 <= j <= ntrace

Input and return

couple:

2D ndarray (or list, or tuple) of size [ncouple , 2] (ex: [[1,1],[1,2],[1,3]] to compute {1,1}, {1,2}, {1,3} cross-correlations

trace:

(int) number of traces in the 2D das section

base:

(int) starting index for trace number, 0 (default) or 1

function return

status:

0 ok; 1 wrong index number (<0 or >ntrace)

Usage example

>>> #For a section containing 500 traces
>>> #Compute correlation of trace 0 and 1 with all other traces
>>> import a1das
>>> ntrace=500
>>> couple =  [list([0,j]) for j in range(0,ntrace)]
>>> couple += [list([1,j]) for j in range(0,ntrace)]
>>> a1das.xcor.register_couple(couple,ntrace=ntrace)
xcor.reset_couple()

Description

Reset the list of couple of trace registered for xcorr computation

xcor.register_par(proc_key, val=[])

Description

Register the parameters used to pre-process data for correlation computation

Input

proc_key:

(str) string referring to a known process

xcor.list_processing()

Description

print the list of pre-processing available when computing cross-correlation

xcor.trace_index_list()

Description

trace_index_list: return the indices of the traces used to compute correlation

xcor.sort_xcorr(i0, ijx, xcorr=None)

Description

I, J = sort_xcorr(i0, ijx) or xcorr_sorted = sort_xcorr(i0, ijx ,xcorr)

return the xcorr indices I that correspond to correlation couple {i0, ?}, xcorr[I,:] are the cross-correlation that implies trace i0 Among these correlation, those with index J need to be time-flipped If an xcorr[ncor, nlags] array is given as input, the function returns a sorted version of xcorr

Input

src:

index of the “source” trace

ijx:

array of indices obtained from the compute_xcorr()

function return

I:

array of index such that xcorr[I, :] correspond to cross-correlation of couple {src, ?} or {?, src)

J:

array of index xcorr[J,:] where one needs to flip time to get {i0,?} xcorr instead of {?,i0}

xcorr_sorted[ncor, ntime]:

sub-array of xcorr that corresponds only to {i0,?} cross-correlations

Usage example

>>> #a1 = a <A1Section> instance
>>> xcorr, lags, ijx, ier = xcor.compute_xcorr(a1,lag=lag,stack=True,verbose=0)
>>>
>>> # return all the indices {i1} that correspond to correlations {i0,?}
>>> I, J = xcor.sort_xcorr(i0,ijx)
>>> xcorrs = xcor.sort_xcorr(i0,ijx, xcorr)
>>> xcorr_sorted = xcorr[I,:]
>>> xcorr_flipped = numpy.flip(xcorr_sorted[J],axis=1)
>>> # xcorrs and xcorr_flipped are identical
xcor.get_nxcorr()

get the number of xcorr couple of trace that have been registered for the computation of cross correlation in the current section