I'm newbie in python and numpy. I have to calculate the meyer wavelet in time domain using spectral analysis. Due to https://academic.oup.com/gji/article/116/1/119/635254 the wavelet can be estimated in frequency domain by
I have tried to write a code in python to construct this wavelet in time domain
import numpy as np
import matplotlib.pyplot as plt
N = 1024
fs = 50
df = fs / N
f = np.arange(-N / 2, N / 2) * df
w = 2 * np.pi * f
def hw(w):
h = []
for i in np.arange(len(w)):
if w[i] > 0:
h.append(np.exp(-1 / w[i] / w[i]))
else:
h.append(0)
return np.array(h)
def gw(w):
return hw(4 * np.pi / 3 - w) / (hw(w - 2 * np.pi / 3) + hw(4 * np.pi / 3 - w))
def phiw(w):
return np.sqrt(gw(w) * gw(-w))
def syw(w):
return np.exp(-1j * w / 2) * np.sqrt(phiw(w / 2) * phiw(w / 2) - phiw(w) * phiw(w))
s = syw(w)
st = np.fft.ifft(np.fft.ifftshift(s))
st = np.fft.fftshift(np.real(st))
t = np.arange(-N / 2, N / 2) / fs
plt.figure(1)
plt.plot(t, st)
plt.xlim([-10, 10])
plt.show()
but I got wrong amplitude in time domain!!??
ps: The domain and amplitude can't be retrived through this code. Any suggestion is welcome
Replace this line:
st = np.fft.fftshift(np.real(st))
with
st = np.fft.fftshift(np.real(st)) * fs
and you will get the correct scaling, I believe.
Dr. Saturn