function [f0slope,sfm_amp,sfm_frq,sfm_phase]=fit_f0(f0data, ploton);
% this function tries to fit f0 as a function of time with a combination of
% a linear fm plus a sinusoidal fm;
% input is an x by 2 matrix with the time values in the first column and
% the corresponding f0 values in the 2nd column
% the matrix must have at least 5 rows
% outputs are four values
% 1st output: linear fm slope (in kHz/ms)
% 2nd output: amplitude of sfm (in kHz)
% 3rd output: frequency of SFM component (in Hz)
% 4th output: phase of SFM component (in rad)

if size(f0data,1)>5,
    % fit a linear regression
    [p,s]=polyfit(f0data(:,1),f0data(:,2),1);
    f0slope=p(1)/1e6;   % in kHz per millisecond
    linfit=polyval(p,f0data(:,1));
    if ploton==1,
        figure(1),
        subplot(2,1,1);
        plot(f0data(:,1),f0data(:,2),f0data(:,1),linfit);
        title(['f0 slope = ' num2str(f0slope) ' kHz/ms']);
    end
    
    % fit a sinusoid to the residual
    x=f0data(:,1)-min(f0data(:,1)); % reset time to start at 0
    y=f0data(:,2)-linfit;
    
    % estimate for sfm_amp
    b0(1)=(max(y)-min(y))/2;
    % estimate for sfm_freq
    b0(2)= 1/(mean(diff(find(diff(sign(y)))))*2*mean(diff(x)));
    % estimate for sfm_phase
    if diff(y(1:2))>0
        b0(3)=0;
    else
        b0(3)=1;
    end
    
    b=nlinfit(x,y,@(b,x) b(1).*sin(2*pi*x*b(2)+b(3)*pi),b0);
    sinfit=b(1).*sin(2*pi*x*b(2)+b(3)*pi);
    sfm_amp=b(1)/1e3;   % SFM Amplitude (kHz)
    sfm_frq=b(2);
    sfm_phase=b(3);
    
    if ploton==1,
        figure(1),
        subplot(2,1,2);
        plot(x,y,'b',  x,sinfit, 'r')
        title(['SFM FRQ = ' num2str(sfm_frq) ' Hz; SFM AMP = ' num2str(sfm_amp) ' kHz']);
    end
else
    f0slope=nan;
    sfm_amp=nan;
    sfm_frq=nan;
    sfm_phase=nan;
end