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:

enter image description here

However, after applying the deconvolution I get something like this:

enter image description here

I'm very confused by those two high spikes at the beginning 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')

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.

  • $\begingroup$ Can I please ask if this was resolved? $\endgroup$ – A_A Jul 7 '20 at 16:29

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.

  • $\begingroup$ Thank you very much for your answer, I've been trying to understand (and apply) what you told me to my code, I'm sorry if I make some obvious mistakes but I'm barely starting into this topic. For Gauss, I defined it with frecg, so: frecg = np.linspace(-numpoints//2, numpoints//2, numpoints); defined the fft using rfft (and without the numpoints//2 range); and used a gauss parameter of 3. The results are still odd. Just for testing, I set the gauss to 1 in erefrec (to see what happens if there is no filter) and I get similar results, with two big spikes at the ends. $\endgroup$ – Danilo Chamorro Riascos Jun 30 '20 at 16:56
  • $\begingroup$ @DaniloChamorroRiascos Thank you for letting me know. No need to apologise. It would be difficult to try and "debug" the situation remotely. My message to you is that you need to understand how the discrete time and frequency work in the discrete Fourier transform space, especially when it comes to defining filters. Beyond this, I could give a closer look provided that I can replicate this problem locally. For this, it would be useful to have a link towards the data. Is this possible? $\endgroup$ – A_A Jul 1 '20 at 11:41

Your Answer

