% E4896_L06_lpc_diary.m
% Diary to show LPC analysis and synthesis
% 2011-02-21, 2012-02-20 Dan Ellis dpwe@ee.columbia.edu

% Read in waveform
[d,sr] = wavread('sm1_cln.wav');
% Resample to 8 kHz, so LPC poles are "focused" on critical sub-4kHz region
d = resample(d,1,2);
sr = sr/2;
% Find a nice piece of vowel waveform
subplot(311)
plot(d)
% .. around sample 1.5 x 10^4.  Pull 256 samples:
dd = d(1.5e4+[1:256]);
plot(dd)
% apply a hanning window
wdd = hanning(256).*dd;
subplot(312)
plot(wdd)
% check its spectrum (note harmonics + broad resonances)
subplot(313)
plotspec(wdd,sr)
% Find LPC: first calculate autocorrelation
rxx = xcorr(wdd);
% let's fit a 12th order LPC
p = 12;
% Just need the p+1 values at 0
rxx = rxx(length(wdd)+[0:p]);
% solve the normal equation 
R = toeplitz(rxx(1:p));
an = inv(R)*rxx(2:(p+1));
% convert into an actual filter
aa = [1 -an'];
% find this filter's spectrum
[H,W] = freqz(1,aa);
% plot alongside the actual spectrum - captures broad shape
hold on; plot(W*sr/2/pi,dB(H)-40,'-r'); hold off
% apply inverse of this filter to original signal to get residual
rdd = filter(aa,1,wdd);
% plot alongside original - note how glottal pulses are sharpened
subplot(312)
hold on; plot(rdd-0.2,'-r'); hold off
% plot the poles on the zplane, note how they match resonant peaks
subplot(311)
zplane(1,aa)

% Short-time LPC analysis-synthesis
% 8th-order LPC modeling
[a,g,e] = lpcfit(d,8);

% Noise-excited resynthesis
rd = lpcsynth(a,g);
% Listen
soundsc(rd,sr)

% Compare spectrograms
subplot(211)
specgram(d,256,sr);
subplot(212)
specgram(rd,256,sr);
% Clearly the same broad shape, but without pitch structure.

% Residual-excited resynthesis (i.e. using the 'e' argument)
rde = lpcsynth(a,g,e);
soundsc(rde,sr)
% Perceptually identical to original
% What is the excitation signal like?
% It should contain everything *not* in the noise-excited reconstruction
soundsc(e,sr)
% Weird, flat power in time and frequency, still intelligible!
specgram(e,256,sr);
% Interesting that formants still visible, resulting from slight 
% modeling mismatch in the LPC

% Compare original waveform with residual
plot(d)
% Choose a portion straddling fricative-vowel in "football"
subplot(311)
plot(d(2000+[1:1000]))
subplot(312)
plot(rde(2000+[1:1000]))
% Residual-excited resynth is close but not identical due to
% the overlap-adding of adjacent frames
subplot(313)
plot(e(2000+[1:1000]))
% Residual is basically noise during fricatives or silence;
% clearly shows pitch pulses during voicing

% Buzz-hiss encoding of residual e
length(e)
% 27840 samples
% Represent as pitch (or noise flag) every 128 samples 
% (from a 256 point window)
P = lpcBHenc(e);
length(P)
% just 217 values

% Decode back into waveform
re = lpcBHdec(P);
% What does it sound like?
soundsc(re,sr)
% .. alternating buzz and hiss
% Use for resynthesis
rdre = lpcsynth(a,g,re);
soundsc(rdre,sr)
% Somewhat unnatural, but not bad!