From 3a5de3fc7b200caf5ffc44765cb1e4b4228cca3f Mon Sep 17 00:00:00 2001
From: Connor Olding <cloningdonor@gmail.com>
Date: Mon, 26 Oct 2015 04:04:32 -0700
Subject: [PATCH] update 13

---
 lib/__init__.py |  1 +
 lib/cepstrum.py | 43 +++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 44 insertions(+)
 create mode 100644 lib/cepstrum.py

diff --git a/lib/__init__.py b/lib/__init__.py
index 4511c25..ef29685 100644
--- a/lib/__init__.py
+++ b/lib/__init__.py
@@ -9,6 +9,7 @@ from .wav import *
 from .planes import *
 from .fft import *
 from .bs import *
+from .cepstrum import *
 
 import numpy as np
 from matplotlib.pylab import show
diff --git a/lib/cepstrum.py b/lib/cepstrum.py
new file mode 100644
index 0000000..daaf1c7
--- /dev/null
+++ b/lib/cepstrum.py
@@ -0,0 +1,43 @@
+import numpy as np
+
+# fast cepstrum and inverted fast cepstrum
+fcs  = lambda s: np.fft.ifft(np.log(np.fft.fft(s)))
+ifcs = lambda s: np.fft.fft(np.exp(np.fft.ifft(s)))
+
+# magnitude
+mcs = lambda s: (np.abs(np.fft.ifft(np.log(np.abs(np.fft.fft(s))**2)))**2)[:len(s)//2]
+
+def clipdb(s, cutoff=-100):
+    as_ = np.abs(s)
+    mas = np.max(as_)
+    if mas == 0 or cutoff >= 0:
+        return s
+    thresh = mas*10**(cutoff/20)
+    return np.where(as_ < thresh, thresh, s)
+
+def fold(r):
+    # via https://ccrma.stanford.edu/~jos/fp/Matlab_listing_fold_m.html
+    # Fold left wing of vector in "FFT buffer format" onto right wing
+    # J.O. Smith, 1982-2002
+    n = len(r)
+    if n < 3:
+        rw = r
+    elif n % 2 == 1:
+        nt = (n + 1)//2
+        rf = r[1:nt] + conj(r[-1:nt-1:-1])
+        rw = np.hstack((r[0], rf, np.zeros(n-nt)))
+    else:
+        nt = n//2
+        rf = np.hstack((r[1:nt], 0)) + np.conj(r[-1:nt-1:-1])
+        rw = np.hstack((r[0], rf, np.zeros(n-nt-1)))
+    return rw
+
+def minphase(s, os=True):
+    # via https://ccrma.stanford.edu/~jos/fp/Matlab_listing_mps_m.html
+    if os:
+        s = np.hstack((s, np.zeros(len(s))))
+    cepstrum = np.fft.ifft(np.log(clipdb(np.fft.fft(s), -100)))
+    signal = np.real(np.fft.ifft(np.exp(np.fft.fft(fold(cepstrum)))))
+    if os:
+        signal = signal[:len(signal)//2]
+    return signal