This python program is based on the work of Mao and Mordret, that you can find here : https://github.com/shujuanmao/dt-wavelet
Contact : Higueret Quentin (quentin.higueret@univ-grenoble-alpes.fr) and Aurélien Mordret (aurelien.mordret@univ-grenoble-alpes.fr)
This package contains codes and test data for measuring seismic travel-time shifts in the time-frequency domain using the wavelet cross-spectrum analysis.
All the files must be in the same folder when you launch the program
Python 3 and the PyCWT packages from https://github.com/regeirk/pycwt are needed to run the codes
You can use pip to install this package :
$ pip install pycwt
-
xwt.py
The core function to calculate dt in the time-frequency domain by wavelet cross-spectrum analysis. -
plotting_example.py
An example of using the cwt function on synthetic data. Plots come with one click. -
ori_waveform.npy
new_waveform.npy
Two synthetic waveforms for testing the codes, ori_waveform and new_waveform
The synthetic seismograms are generated using velocity models by a homogeneous background superimposed with random heterogeneities.
The perturbation between the current and reference velocity models is a 0.05% homogeneous dv/v throughout the medium. (If interested, see Section 3.1 in the following reference for more details.) -
time.npy
fs.npy
The time vector and sampling frequency associated with the synthetic waveforms
WXamp, WXspec, WXangle, Wcoh, WXdt, freqs = xwt(trace_ref,trace_current,fs,ns,nt,vpo,freqmin,freqmax,nptsfreq)
-
Input
trace_ref,trace_current : Two vectors, reference and current time series.
fs : Sampling Frequency // Positive scalar, sampling frequency.
ns : NumScalesToSmooth // Positive integer, indicating the length of boxcar window.
nt : DegTimeToSmooth // Positive scalar, indicating the length of the Gaussian window.
vpo : VoicesPerOctave // Even integer from 4 to 48, indicates how fine the frequency is discretized. Should be no less than 10.
freqmin : The starting value of the frequency vector, in Hz.
freqmax : The ending value of the frequency vector, in Hz.
nptsfreq : Number of frequency samples to generate between the starting and ending value. -
Output
WXamp : Matrix of amplitude product of two CWT in time-frequency domain.
WXspec : Complex-valued matrix, the wavelet cross-spectrum.
WXangle : Matrix of the angle of the complex argument in WXspec.
Wcoh: Matrix of wavelet coherence.
WXdt : Matrix of time difference and phase difference, respectively between the two input time series in time-frequency domain.
freqs : Vector of frequencies used in CWT, in Hz
Using the cwt function and the synthetic data, we can perform a cross-wavelet transform.
This is the plot we obtain using the plotting_example.py program :
(If interested, compare the image obtain here with the Figure 3 in this article Geophysical Journal International, Volume 221, Issue 1, April 2020, Pages 550–568, https://doi.org/10.1093/gji/ggz495)
Implementation of the cone of influence for each frequency bands
S.Mao, A.Mordret, M.Campillo, H.Fang, R.D.van der Hilst, (2019), On the Measurement of Seismic Travel-Time Changes in the Time-Frequency Domain with Wavelet Cross-Spectrum Analysis, GJI, In Review.
Torrence, C. and Compo, G. P.. A Practical Guide to Wavelet Analysis.