Center value of source time value from mt_deconvolve is always zero

[VioletaSeo]

Hi. I am using mtspec library to compute multitaper spectra/spectral ratio in frequency domain. I am also trying to use the function mt_deconvolve to obtain reliable source time function in the time domain. Fortunately, there is good tutorial Multitaper deconvolution for earthquake source studies in the documentation, so I followed the tutorial to compute source time function. Lowpass filter is applied to STF in the documentation to make its shape better (I assume), but I tested various datasets with/without any filtering and found out that the center value of STF is always close to "zero", which makes the shape of STF undesirable. I felt that I am misunderstanding something serious, so I want to ask here what I am doing wrong.

I guess the similar question was raised in https://github.com/krischer/mtspec/issues/26, but I could not get satisfying remedy in that thread.

I attach the code I used below: I used waveform file main.sac which was used in the tutorial. In my case, I tried deconvolving mainshock waveform by itself, because it should lead to delta function in ideal case. In many papers, I have seen that such test is often conducted to check the minimum resolution that can be obtained from deconvolution processing. I hope my issue is reproducible since I am using shared waveform data.

# Read data as ObsPy trace
main_sacfile="./main.sac"
egf_sacfile="./main.sac" # Same file as mainshock

tr_main=read(main_sacfile)[0]
tr_egf=read(egf_sacfile)[0]

# Below, I follow documentation tutorial, but without lowpass filtering
sr=tr_main.stats.sampling_rate
dt=1.0 / sr
tbw=4

r=mt_deconvolve(tr_main.data, tr_egf.data, dt, nfft=tr_main.stats.npts, time_bandwidth=4, number_of_tapers=2*tbw-1, weights="constant", fmax=49, demean=True) # I set fmax=49 here (no cosine-squared filter leads to very noisy result)
decon=r['deconvolved']

# Time vector for RSTF
time=np.arange(0, len(decon))*dt

M=np.arange(0, len(decon))
N=len(M)

SeD=np.where(np.logical_and(M>=0, M<=N/2))
d1=decon[SeD]
SeD2=np.where(np.logical_and(M>N/2, M<=N+1))
d2=decon[SeD2]

# Relative source time function
stf=np.concatenate((d2, d1))

# Plotting
plt.figure(figsize=(16, 4))
plt.plot(time, stf, linewidth=1, color="r")
plt.axhline(y=0.0, color="k", lw=1)
plt.xlabel("Time [sec]", fontsize=25)
plt.ylabel("Amplitude", fontsize=25)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

plt.show();

I attach the result obtained by code above. As it can be clearly seen, there is valley amidst of triangular shaped STF (exactly at the center). I thought such phenomenon should not occur, so I also tried same deconvolution, but after resampling waveform to 1000 Hz (using lanczos interpolation implemented in ObsPy).

# Read data as ObsPy trace
main_sacfile="./main.sac"
egf_sacfile="./main.sac"

tr_main=read(main_sacfile)[0]
tr_egf=read(egf_sacfile)[0]

# Resampling to 1000 Hz
tr_main=tr_main.interpolate(sampling_rate=1000, method="lanczos", a=20)
tr_egf=tr_egf.interpolate(sampling_rate=1000, method="lanczos", a=20)

sr=tr_main.stats.sampling_rate
dt=1.0 / sr

tbw=4

r=mt_deconvolve(tr_main.data, tr_egf.data, dt, nfft=tr_main.stats.npts, time_bandwidth=4, number_of_tapers=2*tbw-1, weights="adaptive", fmax=49, demean=True)
decon=r['deconvolved']
freq=r['frequencies']

# Time vector for RSTF
time=np.arange(0, len(decon))*dt

M=np.arange(0, len(decon))
N=len(M)

SeD=np.where(np.logical_and(M>=0, M<=N/2))
d1=decon[SeD]
SeD2=np.where(np.logical_and(M>N/2, M<=N+1))
d2=decon[SeD2]

# Relative source time function
stf=np.concatenate((d2, d1))

# Plotting
plt.figure(figsize=(16, 4))
plt.plot(time, stf, linewidth=1, color="r")
plt.axhline(y=0.0, color="k", lw=1)
plt.xlabel("Time [sec]", fontsize=25)
plt.ylabel("Amplitude", fontsize=25)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

plt.show();

Interestingly, the valley became much narrower relative to total duration of STF, although it still occupies two samples (going down and up), which means that problem still exists.

To show that in both cases, the durations of STFs exactly coincide regardless of "width of valley", I attach another figure below. I tried three sampling rates (100 Hz (default), 1000 Hz, 10000Hz).

I saw in many papers that this mtspec package is used to calculate STFs in time domain and I have never seen them looking like what I calculated. I also tested water-level deconvolution method and found out that the shape is same, but with no "valley". I would like to really know what I am doing wrong. Perhaps I am misunderstanding something about inverse Fourier transform, but I simply followed the tutorial script.

I tried to also use fortran mtspec code, but unfortunately I am suffering from compiling issue.

Thank you for reading this long text!

++ Actually, I thought the center value was exactly zero at first, but by checking again, it is not exactly zero, but very close to zero. Obviously, I think is should not occur for the simplest case I demonstrated (ex. delta function test).

[VioletaSeo]

I used Python library recently developed by Dr. Prieto (https://github.com/gaprieto/multitaper), and found out that these problems do not occur. I recommend using modules MTSpec and MTCross of multitaper package instead of mt_deconvolve of mtspec for calculating time domain STFs unless you know how to cope with it.

[krischer]

@VioletaSeo Cool - thanks for letting me know about this package! It makes a lot of sense to recommend this new package as it is Python implementation by the author of the fortran module this package here wraps.

The central readme here now recommends usage of this package.