% ex102.m
% get a waveform
fn=input('Enter WAV filename : ','s');
[x,fs,nb]=wavread(fn);
ms2=floor(fs*0.002);
ms10=floor(fs*0.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');
% process in chunks
pos=1;
fx=[];
while (pos+ms30) <= length(x)
    y=x(pos:pos+ms30-1);
    y=y-mean(y);
    % calculate autocorrelation
    r=xcorr(y,ms20,'coeff');    
    r=r(ms20+1:2*ms20+1);
    % search for maximum  between 2ms (=500Hz) and 20ms (=50Hz)
    [rmax,fxval]=max(abs(r(ms2:ms20)));
    if (rmax > 0.5)
        fx=[fx fs/(ms2+fxval-1)];
    else
        fx=[fx 0];
    end;
    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)');