-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathxwt.py
170 lines (124 loc) · 6.67 KB
/
xwt.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
###########################################################################
## XWT Function:
# Using wavelet cross-spectrum to calculate the lapse-time and frequency-dependent
# travel-time changes between two time series.
## USAGE
# WXamp, WXspec, WXangle, Wcoh, WXdt, freqs, coi = 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 to 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.
# !!!! This WXdt is obtained with wrapped phase difference.
# If needed, the user can also produce time difference with
# unwrapped phase from WXangle. !!!!
# freqs : Vector of frequencies used in CWT, in Hz
# coi : Cone of influce. Vector of the maximum period of useful information at each arrival time.
# Note that here coi means MAX PERIODS, in Matlab codes it refers to MIN FREQUENCY.
## EXAMPLE:
# WXamp, WXspec, WXangle, Wcoh, WXdt, freqs, coi = xwt(trace_ref,trace_current,fs,ns,nt,vpo,freqmin,freqmax,nptsfreq)
## Authors: Higueret Quentin ([email protected]) & Aurélien Mordret ([email protected])
# Created: Feb., 2021
## Reference:
# 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
###########################################################################
## Import modules
# The PyCWT package must be installed on the system
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import pyplot, transforms
from scipy.io import loadmat
import pycwt as wavelet
from pycwt.helpers import find
from scipy.signal import convolve2d
import warnings
## Disable Warnings
warnings.filterwarnings('ignore')
## conv2 function
# Returns the two-dimensional convolution of matrices x and y
def conv2(x, y, mode='same'):
return np.rot90(convolve2d(np.rot90(x, 2), np.rot90(y, 2), mode=mode), 2)
## nextpow2 function
# Returns the exponents p for the smallest powers of two that satisfy the relation : 2**p >= abs(x)
def nextpow2(x):
res = np.ceil(np.log2(x))
return res.astype('int')
## Smoothing function
# Smooth the dataset
def smoothCFS(cfs, scales, dt, ns, nt):
N = cfs.shape[1]
npad = 2 ** nextpow2(N)
omega = np.arange(1, np.fix(npad / 2) + 1, 1).tolist()
omega = np.array(omega) * ((2 * np.pi) / npad)
omega_save = -omega[int(np.fix((npad - 1) / 2)) - 1:0:-1]
omega_2 = np.concatenate((0., omega), axis=None)
omega_2 = np.concatenate((omega_2, omega_save), axis=None)
omega = np.concatenate((omega_2, -omega[0]), axis=None)
# Normalize scales by DT because we are not including DT in the angular frequencies here.
# The smoothing is done by multiplication in the Fourier domain.
normscales = scales / dt
for kk in range(0, cfs.shape[0]):
F = np.exp(-nt * (normscales[kk] ** 2) * omega ** 2)
smooth = np.fft.ifft(F * np.fft.fft(cfs[kk - 1], npad))
cfs[kk - 1] = smooth[0:N]
# Convolve the coefficients with a moving average smoothing filter across scales.
H = 1 / ns * np.ones((ns, 1))
cfs = conv2(cfs, H)
return cfs
## xwt function
def xwt(trace_ref, trace_current, fs, ns, nt, vpo, freqmin, freqmax, nptsfreq):
# Choosing a Morlet wavelet with a central frequency w0 = 6
mother = wavelet.Morlet(6.)
# nx represent the number of element in the trace_current array
nx = np.size(trace_current)
x_reference = np.transpose(trace_ref)
x_current = np.transpose(trace_current)
# Sampling interval
dt = 1 / fs
# Spacing between discrete scales, the default value is 1/12
dj = 1 / vpo
# Number of scales less one, -1 refers to the default value which is J = (log2(N * dt / so)) / dj.
J = -1
# Smallest scale of the wavelet, default value is 2*dt
s0 = 2 * dt # Smallest scale of the wavelet, default value is 2*dt
# Creation of the frequency vector that we will use in the continuous wavelet transform
freqlim = np.linspace(freqmax, freqmin, num=nptsfreq, endpoint=True, retstep=False, dtype=None, axis=0)
# Calculation of the two wavelet transform independently
# scales are calculated using the wavelet Fourier wavelength
# fft : Normalized fast Fourier transform of the input trace
# fftfreqs : Fourier frequencies for the calculated FFT spectrum.
cwt_reference, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(x_reference, dt, dj, s0, J, mother, freqs=freqlim)
cwt_current, _, _, _, _, _ = wavelet.cwt(x_current, dt, dj, s0, J, mother, freqs=freqlim)
scales = np.array([[kk] for kk in scales])
invscales = np.kron(np.ones((1, nx)), 1 / scales)
cfs1 = smoothCFS(invscales * abs(cwt_reference) ** 2, scales, dt, ns, nt)
cfs2 = smoothCFS(invscales * abs(cwt_current) ** 2, scales, dt, ns, nt)
crossCFS = cwt_reference * np.conj(cwt_current)
WXamp = abs(crossCFS)
# cross-wavelet transform operation with smoothing
crossCFS = smoothCFS(invscales * crossCFS, scales, dt, ns, nt)
WXspec = crossCFS / (np.sqrt(cfs1) * np.sqrt(cfs2))
WXangle = np.angle(WXspec)
Wcoh = abs(crossCFS) ** 2 / (cfs1 * cfs2)
pp = 2 * np.pi * freqs
pp2 = np.array([[kk] for kk in pp])
WXdt = WXangle / np.kron(np.ones((1, nx)), pp2)
return WXamp, WXspec, WXangle, Wcoh, WXdt, freqs, coi