pythonnumpyfftwavelet

Some problem in calculating the meyer wavelet using IFFT


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 enter image description here

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!!??

enter image description here

enter image description here

ps: The domain and amplitude can't be retrived through this code. Any suggestion is welcome

The figure in the article


Solution

  • 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