cepst.m

cepstral analisys for pitch and F0 detection

function fx = cepst (x,fs)
    x = x';
    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=ifft(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)');
end