I'm working on some code for receiver function method in seismology. For anyone one not into the topic, it's just a deconvolution of two time series (seismograms). This can be done in the time domain as well as in the frequency domain, however, the method I'm working on is doing it in the latter. The function I'm applying is this one, as described in Langston, 1979:
The phi term is what we call water level, and what it does is that it "fills" the low frequencies in the DFT so it can ignore the noise and smooth the result. The gaussian filter is also intended for smoothing of the result.
Both R and V seismograms looks like this, an usual seismogram:
However, after applying the deconvolution I get something like this:
I'm very confused by those two high spikes at the begining and the end of the result, it should be smooth and not spiky.
This is my code:
verfrec = sp.fft(stv.data, numpoints)[:numpoints//2]
radfrec = sp.fft(stradial.data, numpoints)[:numpoints//2]
frec = np.linspace(1, numpoints//2, numpoints//2)
gauss = np.exp((-frec**2)/((4*60)**2))
waterlevel = np.abs(verfrec * np.conjugate(verfrec))
waterlevel = np.maximum(waterlevel, 0.1*np.max(waterlevel))
erefrec = radfrec * gauss * np.conjugate(verfrec) / waterlevel
inverse = sp.ifft(erefrec)
fig = plt.figure(figsize=(12, 4))
ax = fig.add_subplot(1, 1, 1)
ax.plot(frec, inverse, 'sandybrown')
plt.title('Sismograma')
fig.show()
I'm splitting the transforms range in the middle since the real input will produce symmetric outputs. The 60 in the gauss definition is the gauss parameter and it should be between 0.1 and 1.5 or so, with those values I get wrong results so I'm guessing its a magnitude thing.
Any help is very appreciated.