Why is my Cepstrum mirrored and how do I relate it back to the original signal to find an echo

Question

I am not a signal processor but need to use Cepstrums in a bigger project I am writing to help identify echoes in seismic signals.

I need to understand why the Real and Complex Cepstrums I create, in Python, have the right hand side of the output being the mirror image of the left hand side. I get a similar result when I create the Complex Cepstrum - although here the mirroring effect is both left to right and up to down. I have based my code on numerous published articles available in the public domain (e.g the paper by R.B. Randall).

I have tried using both synthetic data and real seismic data. But all the results have this mirroring. For simplicity I show the results derived from the synthetic signal.

The first two images show the input synthetic data and its associated spectra.

The next image show the Real Cepstrum. Here the mirroring effect is clear, with the right hand side of the output being the mirror image of the left hand side.

The last two images show the result of the (complex) frequency spectra and then the complex Cepstrum. Here the mirroring effect is a combination of flipping the right hand side both right to left and up to down.

What am I doing wrong? The mirroring is baffling me. Also I know the Cepstrum should be the same length as the input signal but taking the first half of the results I get means the Cepstrum is too short.

Once I have the Cepstrum what is the best way to then us sit to determine the echo in the original signal from the Cepstrum ?

Below is the Scipy / Python code I am using to derive the results shown above

el3 is the input signal

SampRate =20

---

REAL CEPSTRUM

yfft = rfft(el3)
xfft = rfftfreq(len(el3), 1 / SampRate)
yfft_abs = (np.abs(yfft))
CEPsig = irfft(np.log(yfft_abs))

---

COMPLEX CEPSTRUM

A_f = fft(el3)
A_fft = fftfreq(len(el3), 1 / SampRate)
C_CEPsig = ifft(np.log(A_f))

---

Code for plotting cepstrum is

plt.figure(figsize=(30,4))
plt.ylim(-1,1)
plt.xticks(np.arange(0, len(CEPsig), 100))
plt.plot(CEPsig, color = "red")
plt.title("Real Cepstrum of synthetic data", size = 20)
plt.show()

Code for calculating magnitude and phase (in radians) from Complex Cepstrum is

C_CEPsig_abs = np.abs(C_CEPsig)
C_CEPsig_ang = np.angle(C_CEPsig)

Code for plotting magnitude and phase is

plt.figure(figsize=(30,4))
plt.plot( C_CEPsig_ang, color = "black")
plt.title("Phase of Cepstrum", size = 20)
plt.show()

plt.figure(figsize=(30,4))
plt.plot( C_CEPsig_abs, color = "red")
plt.title("Magnitude of Cepstrum", size = 20)
plt.show()

These plots are shown below. Although I think the Phase vs Quefrency plot may be OK (?) I still see the mirroring effect on the Magnitude of the Complex Cepstrum vs the Quefrency plot.

As an aside from the code below the number of samples in the quefrency_vector_impulse (= 901) is half that of the number of samples in the C_CEPsig (= 1800)

A_fft = fftfreq(len(el3), 1 / SampRate)
df = A_fft[1] - A_fft[0]
quefrency_vector_impulse = rfftfreq(C_CEPsig.size, df)

Answer 1

Cepstrum is a Fourier transform. Since \begin{align} \sin(\omega_i t)&=-\sin(-\omega_i t)\\\cos(\omega_i t)&=\cos(-\omega_i t) \end{align} the coefficients representing the magnitude of those components are going to reflect this behavior. e.g symmetric for the cosine and anti-symmetric for the sine.

Answer 2

Let $\log$ denote the natural logarithm, $\mathcal{F}$ denote the Fourier transform, and $\mathcal{F}^{-1}$ denote the inverse Fourier transform. You can also use a different logarithm base (i.e. 10), depending on the scaling of your data. There is no single "correct" way to select the base for every situation. The definition of the real cepstrum is $$C_r=\mathcal{F}^{-1}\Big(\log\left(\big|\mathcal{F}\big(f(t)\big)\big|\Big)\right).$$

In contrast, the complex cepstrum is defined as $$C_c=\mathcal{F}^{-1}\Big(\log\left(\mathcal{F}\big(f(t)\big)\Big)\right).$$

Note that these two differ slightly from the more often used power cepstrum $$C_p=\Big|\mathcal{F}^{-1}\Big(\log\left(\big|\mathcal{F}\big(f(t)\big)\big|^2\Big)\right)\Big|^2.$$

As Peter mentioned, Fourier transforms (and their inverses) of real inputs have a complex conjugate symmetry between positive and negative frequencies. Currently, you are directly plotting the complex values, which in most plotting packages outright discards the imaginary part. Hence the symmetry in your plot.

I encourage you to plot the magnitude and phase (numpy.abs and numpy.angle, respectively) as a function of quefrency of the complex quantities calculated above. Recall the following relations for complex numbers: $$z=a+ib$$ $$|z|=\sqrt{a^2+b^2}$$ $$\angle z = \tan^{-1}\left(\frac{b}{a}\right)$$ where the arctan domain is kept unrestricted through the use of atan2 in most programming languages.

Answer 3

The solution I have found to ensure that the mirroring (symmetry) is not observed in a Real Cepstrum in Python is to double the length of the input. Python allows you to do this with the variable 'n' as in the following code line

yfft = rfft(el4, n=len(el4)*2 )

The will pad the input with zeros. So when the Spectra and Cepstrum are created the number of their samples is double. The second half of the Cepstrum is mirrored but will not overlap with the first half. Taking the first half of the Cepstrum then gives the response required