%pylab --no-import-all inline
import matplotlib.pyplot as plt

from scikits.talkbox import lpc, segment_axis, lpcres, slfilter
from scipy.io import wavfile
from scipy.signal import lfilter, cwt, ricker
from pylab import specgram
from pywt import wavedec, Wavelet

# Load a test speech file (this one was downloaded from http://freesound.org/people/kenders2000/sounds/175234/)
(fs, s) = wavfile.read('test16k.wav')

plt.plot(s)

# Plot 200 time samples
plt.stem(s[16000:16250])

framelen=int(20e-3*fs)
specgram(s,NFFT=framelen, Fs=fs, noverlap=framelen/2)

# Segment the input vector into non-overlapping frames of 20 ms
X = segment_axis(s, framelen)

# Compute the 10th order LPC coefficients (a), the total residual (e), and the reflection coefficients (k) for each frame
(a, e, k) = lpc(X, 20)
# Compute the LPC residual vectors
res = lpcres(X,20)
# Resynthesize first frame
y = lfilter([1.0], a[100], res[100])

plt.subplot(211)
plt.plot(X[100],'b',y,'rx')
plt.subplot(212)
plt.plot((y - X[100])**2)

# Resynthesize the whole file
Y = slfilter(np.ones((len(a),1)),a,X)
yall = np.concatenate(Y)

from audiodisplay import Audio
Audio(data=s, rate=fs)

Audio(data=(yall/max(abs(yall))*2**15).astype(np.int16), rate=fs)

plt.plot(np.concatenate(res))

# Plotting some Daubechies wavelet functions with different widths
wavelet = Wavelet('db2')

for i in range(5):
    phi, psi, x = wavelet.wavefun(level=i+1)
    plt.subplot(5,2,2*i)
    plt.plot(phi)
    plt.subplot(5,2,2*i+1)
    plt.plot(psi)

# 3 level decomposition of the speech signal
coeffs = wavedec(s, 'db2', level=3)

for i in range(len(coeffs)):
    plt.subplot(len(coeffs),1,i)
    plt.plot(coeffs[i])

Audio(coeffs[3].astype(np.int16), rate=fs/2)

Audio(coeffs[2].astype(np.int16), rate=fs/4)

Audio(coeffs[1].astype(np.int16), rate=fs/8)

Audio(coeffs[0].astype(np.int16), rate=fs/8)