Cepstrum

Having a quick look a the cepstrum.

Reproducing this paper of mine.

Prelims



In [1]:

    
import numpy as np
from scipy.fftpack import fft, rfft, irfft, fftfreq, rfftfreq
import scipy.signal
import matplotlib.pyplot as plt
%matplotlib inline



In [2]:

    
# Create RC series
dt = 0.001    # s, equiv. 1 ms
T = 1.0       # s
thick = 0.040 # s, equiv. 40 ms
rc = 0.3

t = np.linspace(0, T, T/dt)
rcs = np.zeros_like(t)
rcs[1000*(T-thick)/2] = rc
rcs[-1000*(T-thick)/2] = -rc

ax1 = plt.subplot(111)
ax1.plot(t, rcs)
ax1.set_ylim(-0.4, 0.4)
plt.show()



In [3]:

    
RCS = fft(rcs)
freq = fftfreq(rcs.size, d=dt)

keep = freq>=0 # only positive frequencies
RCS = RCS[keep]
freq = freq[keep]



In [4]:

    
plt.figure(figsize=(15,5))
ax1 = plt.subplot(111)
ax1.plot(freq, np.abs(RCS))
ax1.set_xlim(0,200)
plt.show()

Basic signal and plot



In [5]:

    
def ricker(f, duration=0.512, dt=0.001):
    t = np.linspace(-duration/2, (duration-dt)/2, duration/dt)
    y = (1.0 - 2.0*(np.pi**2)*(f**2)*(t**2)) * np.exp(-(np.pi**2)*(f**2)*(t**2))
    return t, y



In [6]:

    
dt = 0.001    # s, equiv. 1 ms
T = 0.256     # s
thick = 0.023 # s, equiv. 23 ms
rc = 0.3

t = np.linspace(0, T, T/dt)
rcs = np.zeros_like(t)
rcs[1000*(T-thick)/2] = rc
rcs[-1000*(T-thick)/2] = -rc

tw, w = ricker(f=43., length=0.064, dt=dt)
s = np.convolve(rcs, w, mode='same')

ax1 = plt.subplot(111)
ax1.plot(t, s)
ax1.set_ylim(-0.4, 0.4)
plt.show()



In [7]:

    
plt.plot(tw,w)









    Out[7]:





[<matplotlib.lines.Line2D at 0x10e4c0850>]



In [8]:

    
S = fft(s)
freq = fftfreq(s.size, d=dt)

keep = freq>=0 # only positive frequencies
Sk = S[keep]
freqk = freq[keep]



In [9]:

    
plt.figure(figsize=(15,5))
ax1 = plt.subplot(111)
ax1.plot(freqk, np.abs(Sk))
ax1.set_xlim(0, 250)
plt.show()



In [10]:

    
freqk = freqk[:]
Sk = np.abs(Sk)



In [11]:

    
Q = fft(Sk)
quef = fftfreq(Sk.size, d=freqk[1]-freqk[0])

keep = quef>=0 # only positive frequencies
Qk = Q[keep]
quefk = quef[keep]



In [12]:

    
quefk









    Out[12]:





array([ 0.   ,  0.002,  0.004,  0.006,  0.008,  0.01 ,  0.012,  0.014,
        0.016,  0.018,  0.02 ,  0.022,  0.024,  0.026,  0.028,  0.03 ,
        0.032,  0.034,  0.036,  0.038,  0.04 ,  0.042,  0.044,  0.046,
        0.048,  0.05 ,  0.052,  0.054,  0.056,  0.058,  0.06 ,  0.062,
        0.064,  0.066,  0.068,  0.07 ,  0.072,  0.074,  0.076,  0.078,
        0.08 ,  0.082,  0.084,  0.086,  0.088,  0.09 ,  0.092,  0.094,
        0.096,  0.098,  0.1  ,  0.102,  0.104,  0.106,  0.108,  0.11 ,
        0.112,  0.114,  0.116,  0.118,  0.12 ,  0.122,  0.124,  0.126])



In [13]:

    
Qk[:15]









    Out[13]:





array([ 49.02079808 +0.j        ,  38.79322071-26.98139188j,
        14.95831306-39.5664094j ,  -7.97162925-34.15940996j,
       -19.34154651-18.49851932j, -18.10010185 -3.46717219j,
        -9.86946890 +4.37154908j,  -1.26472664 +4.19896514j,
         3.14791396 -1.28008305j,   1.43459047 -7.77417801j,
        -5.30839023-10.69441329j, -12.90552489 -7.27058113j,
       -16.20291978 +1.23488393j, -12.78345190 +9.9574107j ,
        -4.84207512+14.07714464j])



In [14]:

    
plt.figure(figsize=(15,5))
ax1 = plt.subplot(111)
ax1.plot(quefk, np.abs(Qk))
plt.show()

The big diagram



In [29]:

    
import agilegeo.wavelet as ag



In [76]:

    
dt = 0.001    # s, equiv. 1 ms
T = 0.256     # s
thick = 0.021 # s, equiv. 21 ms
rc = 0.3

# RC SERIES
t = np.linspace(0, T, T/dt)
rcs = np.zeros_like(t)
rcs[1000*(T-3*thick)/2] = rc
rcs[-1000*(T-3*thick)/2] = -rc
rcs[1000*(T-2*thick)/2] = -rc
rcs[-1000*(T-2*thick)/2] = rc
rcs[1000*(T-thick)/2] = rc
rcs[-1000*(T-thick)/2] = -rc

# WAVELET
duration = 0.256
f = f=[10,16,64,80]
#tw, w = ricker(f=32., duration=0.128, dt=dt)
w = ag.ormsby(f=f, duration=duration, dt=dt)
tw = np.linspace(-duration/2, (duration-dt)/2, duration/dt)

# CONVOLVE
s = np.convolve(rcs, w, mode='same')

# Plot it
plt.figure(figsize=(20,3))

ax1 = plt.subplot(131)
ax1.plot(tw, w)
ax1.set_ylim(-0.5, 1.1)

ax1 = plt.subplot(132)
ax1.plot(t, rcs)
ax1.set_xlim(0, 0.256)
ax1.set_ylim(-0.4, 0.4)

ax1 = plt.subplot(133)
ax1.plot(t, s)
ax1.set_xlim(0, 0.256)
ax1.set_ylim(-0.4, 0.4)

plt.show()



In [70]:

    
W = fft(w)
freqw = fftfreq(w.size, d=dt)

RCS = fft(rcs)
S = fft(s)
freq = fftfreq(s.size, d=dt)

keep = freq>=0
keepw = freqw>=0

Wk = W[keepw]
freqwk = freqw[keepw]
RCSk = RCS[keep]
Sk = S[keep]
freqk = freq[keep]

plt.figure(figsize=(20,3))

ax1 = plt.subplot(131)
ax1.plot(freqwk, np.abs(Wk))
ax1.set_xlim(0,200)

ax2 = plt.subplot(132)
ax2.plot(freqk, np.abs(RCSk))
ax2.set_xlim(0,200)

ax3 = plt.subplot(133)
ax3.plot(freqk, np.abs(Sk))
ax3.set_xlim(0, 200)

plt.show()



In [71]:

    
Wk = np.abs(Wk)
RCSk = np.abs(RCSk)
Sk = np.abs(Sk)

Qw = fft(Wk)
Qr = fft(RCSk)
Qs = fft(Sk)

quefw = fftfreq(Wk.size, d=freqwk[1]-freqwk[0])
quefr = fftfreq(RCSk.size, d=freqk[1]-freqk[0])
quefs = fftfreq(Sk.size, d=freqk[1]-freqk[0])

keepw = quefw>=0 # only positive frequencies
Qwk = Qw[keepw]
quefwk = quefw[keepw]
keep = quefr>=0 # only positive frequencies
Qrk = Qr[keep]
Qsk = Qs[keep]
quefk = quef[keep]

plt.figure(figsize=(15,3))

ax1 = plt.subplot(131)
ax1.plot(quefwk, np.abs(Qwk))

ax2 = plt.subplot(132)
ax2.plot(quefk, np.abs(Qrk))

ax3 = plt.subplot(133)
ax3.plot(quefk, np.abs(Qsk))

plt.show()

This does not quite look right... I wonder if I went wrong somewhere. It looks much cleaner than my effort in Excel. But the first cepstral peak should be at about 21 ms (given this thickness) — and it's actually appearing at double that, roughly.



In [ ]: