% ex103.m
% get a waveform
fn=input('Enter WAV filename : ','s');
[x,fs,nb]=wavread(fn);
ms10=floor(fs*.01);
ms30=floor(fs*0.03);
ncoeff=2+fs/1000;           % rule of thumb for formant estimation
% plot waveform
t=(0:length(x)-1)/fs;
subplot(2,1,1);
plot(t,x);
legend('Waveform');
xlabel('Time (s)');
ylabel('Amplitude');
% process in chunks of 30ms
pos=1;
fm=[];      % formant peaks
ft=[];      % formant times
while (pos+ms30) <= length(x)
    y=x(pos:pos+ms30-1);
    y=y-mean(y);
    % find LP filter
    a=lpc(y,ncoeff);
    % find roots
    r=roots(a);
    r=r(imag(r)>0.01);
    ffreq=sort(atan2(imag(r),real(r))*fs/(2*pi));
    for i=1:length(ffreq)
        fm = [fm ffreq(i)];
        ft = [ft pos/fs];
    end
    pos=pos+ms10;
end;
% plot formants trace
subplot(2,1,2);
plot(ft,fm,'+');
legend('Formants');
xlabel('Time (s)');
ylabel('Frequency (Hz)');