123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657 |
- 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
|