using PyCall
matplotlib = pyimport("matplotlib")
PyDict(matplotlib["rcParams"])["figure.figsize"] = (12, 5)
using PyPlot

import SPTK
using WAV
using DSP
using MelGeneralizedCepstrums

x, fs = wavread(joinpath(Pkg.dir("SPTK"), "examples", "test16k.wav"), format="native")
x = convert(Vector{Float64}, vec(x))

# Visualize the speech signal in time-domain
plot(1:endof(x), x, label="a speech signal")
xlim(1, endof(x))
xlabel("sample")
legend()

# Pick a short segment
pos = 3000
fftlen = 1024

# Note that mel-generalized cepstrum analysis assumes power-normalized window so that Σₙ w(n)² = 1
win = DSP.blackman(fftlen) ./ sqrt(sumabs2(DSP.blackman(fftlen))) # or you can write simply `SPTK.blackman(fftlen)`
@assert isapprox(sumabs2(win), 1.0)

xw = x[pos+1:pos+fftlen] .* win

plot(1:endof(xw), xw, linewidth="2", label="a windowed speech signal")
xlim(1, endof(xw))
xlabel("sample")
legend()

# Plotting utility for visualizing spectral envelope estimate
function pplot(sp, envelope; title="envelope")
    plot(sp, "b-", linewidth="2", label="Original log spectrum 20log|X(ω)|")
    plot(20/log(10)*(envelope), "r-", linewidth="3", label=title)
    xlim(1, length(sp))
    xlabel("frequency bin")
    ylabel("log amplitude")
    legend()
end

# Compute spectrum 20log|X(ω)| for a windowed signal
sp = 20log10(abs(rfft(xw)));

# Linear Cepstrum
c = estimate(LinearCepstrum(20), xw)
pplot(sp, real(mgc2sp(c, fftlen)), title="Linear frequency cepstrum based envelope")

# Mel-Cepstrum
mc = estimate(MelCepstrum(20, 0.41), xw)
pplot(sp, real(mgc2sp(mc, fftlen)), title="Mel-cepstrum based envelope")

# LPC Cepstrum 
mgc = estimate(AllPoleCepstrum(20), xw)
pplot(sp, real(mgc2sp(mgc, fftlen)), title="LPC cepstrum based envelope")

# Warped LPC
mgc = estimate(MelGeneralizedCepstrum(20, 0.41, -1.0), xw)
pplot(sp, real(mgc2sp(mgc, fftlen)), title="Warped LPC based envelope")

# Generalized Cepstrum
mgc = estimate(GeneralizedCepstrum(20, -0.35), xw)
pplot(sp, real(mgc2sp(mgc, fftlen)), title="Generalized cepstrum based envelope")

# Mel-Generalized Cepstrum
mgc = estimate(MelGeneralizedCepstrum(20, 0.41, -0.35), xw)
pplot(sp, real(mgc2sp(mgc, fftlen)), title="Mel-generalized cepstrum based envelope")