% 2001-02-27-cep.diary
% Diary of how the cepstrum figures were produced
% dpwe@ee.columbia.edu

[d,r]=wavread('mpgr1_sx419.wav');
% Extract a little bit of vowel
plot(d(7000+[1:2000]))
dx = hanning(512).*d(8500+[1:512]);
subplot(411)
plot(dx)
X = fft(dx);
ff = [0:255]/256 * 8000;
subplot(412)
plot(ff, abs(X(1:256)));
% Just look at bottom half of spec
ff = [0:127]/256 * 8000;
plot(ff, abs(X(1:128)));
% Log spec
subplot(413)
plot(ff, 20*log10(abs(X(1:128))));
axis([0 4000 -50 10])
% Cepstrum is ifft of log mag spec
c = real(ifft(log(abs(X))));  % Real removes v.small imag partg
subplot(414)
plot(c)
% dominated by c(0)
plot(3:200,c(3:200))

% Lifter selects just first 50 points, but retain start and end
% of spec (so remains a pure-real IFT of a pure-real spectrum)
lf = [ones(1,50), zeros(1,413), ones(1,49)];
plot(3:200,c(3:200),1:200,lf(1:200))
lc = lf.*c';
% Go back to the *linear* spectrum this represents
lX = real(exp(fft(lc)));  % Again, will be real anyway
% Show liftered log spect
subplot(413)
plot(ff, 20*log10(abs(X(1:128))),ff,20*log10(lX(1:128)));
axis([0 4000 -50 10])
% Clearly a smoothed version
% Plot on linear too...
subplot(412)
plot(ff, abs(X(1:128)),ff,lX(1:128));
% Minimum phase IR from cepstral utility:
[xh,yh]=rceps(ifft(lX));
subplot(411)
plot(1:512,dx,1:512,yh)
% More or less plausible.  Show that it has the right spec:
subplot(413)
Y = fft(yh);
plot(ff, 20*log10(abs(Y(1:128))),ff,20*log10(abs(lX(1:128))));
% They match