1
$\begingroup$

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

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

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

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.

enter image description here

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)
$\endgroup$
4
  • $\begingroup$ Can you show us your code for calculating and plotting the Cepstrum? $\endgroup$
    – Ash
    Jul 24, 2022 at 19:34
  • $\begingroup$ Welcome to SE.SP! As Ash says, please give us some (pseudo-)code to show us how you're calculating the cepstrum (real and complex). If the FFT is involved, then that will automatically give conjugate symmetric spectrum for real-valued input which will generate the mirroring you're seeing. $\endgroup$
    – Peter K.
    Jul 24, 2022 at 19:59
  • $\begingroup$ Hi Peter, thanks for your response. I have now added the python code that I have used to create the (real and complex) cepstrums. I am using SciPy for the ffts $\endgroup$
    – Chara
    Jul 25, 2022 at 12:49
  • $\begingroup$ Hi Ash, thanks for your help. I have added the code for generating the cepstums and have now added the code for plotting the cepstrum (Matplotlib). $\endgroup$
    – Chara
    Jul 25, 2022 at 12:55

3 Answers 3

0
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Hi Ash, Thanks for your direction. I have modified the code as per your guidance. When I plot the Magnitude and Phase against Quefrency the phase plot looks reasonable(?) but the magnitude looks very akin to the Real Cepstrum (as expected) I have added code to calculate magnitude and angle. On the magnitude of Complex Cepstrum vs Quefrency I still see the mirroring. So I still do not understand where I am going wrong. where am I coding incorrectly? $\endgroup$
    – Chara
    Jul 27, 2022 at 11:07
0
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ Input length has nothing to do with the mirroring, you solved it because you used rfft. If you're taking half of rfft, you're losing information; this almost certainly isn't desired. $\endgroup$ Aug 7, 2022 at 17:56
  • $\begingroup$ I am not taking half the rfft. In fact I am doubling the length of the original input signal. This additional length is padded with zeros. I then apply the rfft to the new signal ( = original signal plus a padding of the same length as the original length but consisting of zeros ) I continue with this doubled signal length for the rest of the cepstrum. How would you deal with the mirroring/symmetry? $\endgroup$
    – Chara
    Aug 8, 2022 at 15:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.