Receiver function, frecuency domain deconvolution not giving logic results

Question

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.

Answer 1

The "spikes" that you see are coming from a "half-done" inverse Discrete Fourier Transform (DFT).

The main thing to note is that when you specify a filter in the frequency domain (as it is hapenning here), that filter has to "fit" the DFT spectrum ordering too.

You conclude (correctly) "I'm splitting the transforms range in the middle since the real input will produce symmetric outputs. " but that is not what is expressed in the code.

If you use numpy.fft.fft() then you have to specify a symmetric filter. Which is clearly not what gauss will be, because it is defined over freq whose linspace goes from 1 (why one?) to numpoints//2. If you plot this, you will see that it is only half the Gaussian and it has to be extended with its mirror image (since you are using numpy.fft.fft()) to include the negative frequencies part.

What you describe (that you only specify half the spectrum because the rest is symmetric) is what the pair numpy.fft.rfft(), numpy.fft.irfft() implement. But if you choose this way of doing it, you have to get rid of numpoints//2 because it is already taken care for you.

Having done this, we then turn to the "cut off" (or equivalent cut off) of the filter. The formulas in Langston79 are in angular frequency and the "standard deviation" (the width of the Gaussian filter) is specified at some $4 \alpha$, which is something like 6Hz (for an $a=1.5$). This has to be defined appropriately on the discrete frequency axis. If you have (for example) an $Fs=1000Hz$ and you are looking for 6 Hz in a 500 point signal, then after the DFT, each bin is $\frac{Fs}{N_{FFT}}$ Hz wide (or 2 Hz in this case). So your discrete $\alpha$ (let's call it $\alpha_d$) is now $3$. If your signal is very long (which I think it is, here, just looking at the duration of those signals), the $\alpha_d$ will work out much higher of course, because the resolution of the DFT is now higher (to still achieve the desired $4 \alpha$).

Hope this helps.