Having a quick look a the cepstrum.
Reproducing this paper of mine.
import numpy as np
from scipy.fftpack import fft, rfft, irfft, fftfreq, rfftfreq
import scipy.signal
import matplotlib.pyplot as plt
%matplotlib inline
# 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()
RCS = fft(rcs)
freq = fftfreq(rcs.size, d=dt)
keep = freq>=0 # only positive frequencies
RCS = RCS[keep]
freq = freq[keep]
import matplotlib.gridspec as gridspec
gs = gridspec.GridSpec(1, 2,width_ratios=[1,2])
plt.figure(figsize=(15,5))
ax1 = plt.subplot(gs[0])
ax1.plot(t, rcs)
ax1.set_ylim(-0.4, 0.4)
ax2 = plt.subplot(gs[1])
ax2.plot(freq, np.abs(RCS))
ax2.set_xlim(0,200)
plt.savefig("/home/matt/images/cepstrum.png")
plt.show()
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
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()
plt.plot(tw,w)
[<matplotlib.lines.Line2D at 0x10e4c0850>]
S = fft(s)
freq = fftfreq(s.size, d=dt)
keep = freq>=0 # only positive frequencies
Sk = S[keep]
freqk = freq[keep]
plt.figure(figsize=(15,5))
ax1 = plt.subplot(111)
ax1.plot(freqk, np.abs(Sk))
ax1.set_xlim(0, 250)
plt.show()
freqk = freqk[:]
Sk = np.abs(Sk)
Q = fft(Sk)
quef = fftfreq(Sk.size, d=freqk[1]-freqk[0])
keep = quef>=0 # only positive frequencies
Qk = Q[keep]
quefk = quef[keep]
quefk
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])
Qk[:15]
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])
plt.figure(figsize=(15,5))
ax1 = plt.subplot(111)
ax1.plot(quefk, np.abs(Qk))
plt.show()
import agilegeo.wavelet as ag
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()
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()
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.