% ex101.m
% get a waveform
fn=input('Enter WAV filename : ','s');
[x,fs,nb]=wavread(fn);
ms2=floor(fs*0.002);
ms10=floor(fs*.01);
ms20=floor(fs*0.02);
ms30=floor(fs*0.03);
% plot waveform
t=(0:length(x)-1)/fs;
subplot(2,1,1);
plot(t,x);
legend('Waveform');
xlabel('Time (s)');
ylabel('Amplitude');
% get window
w=hamming(ms30);
pos=1;
fx=[];
while (pos+ms30) <= length(x)
    y=x(pos:pos+ms30-1);
    % do fourier transform of windowed signal
    Y=fft(y.*w);
    % cepstrum is DFT of log spectrum
    C=fft(log(abs(Y)+eps));
    % search for maximum  between 2ms (=500Hz) and 20ms (=50Hz)
    [c,fxval]=max(abs(C(ms2:ms20)));
    fx=[fx fs/(ms2+fxval-1)];
    pos=pos+ms10;
end;
% plot FX trace
t=(0:length(fx)-1)*0.01
subplot(2,1,2);
plot(t,fx);
legend('FX Trace');
xlabel('Time (s)');
ylabel('Frequency (Hz)');