tomwalters@0: %
tomwalters@0: %    Calculation of F0 from Summary SAI 
tomwalters@0: %    Irino, T.
tomwalters@0: %    31 Jan 00
tomwalters@0: %    modified from CalPeakSumSAI(SAIval, fs, NumPeak);
tomwalters@0: %
tomwalters@0: %    function [F0estSAI, MeanSAI, ModMeanSAI] ...
tomwalters@0: %             = CalF0estSAI(SAIval, fs, Frs, StrobeLoc, MaxF0, PreF0);
tomwalters@0: %    INPUT:  SAIval : SAI value
tomwalters@0: %            fs     : Sampling frequency
tomwalters@0: %            StrobeLoc: Number of Strobe Location (TimeInterval == 0)
tomwalters@0: %	     MaxF0:  Maximum value of F0
tomwalters@0: %	     PreF0:  Pre-Estimated value of F0
tomwalters@0: %    OUTPUT: F0estSAI: Estimated F0 from SAI
tomwalters@0: %            MeanSAI:   Summary of SAI
tomwalters@0: %            ModMeanSAI: Modified version of MeanSAI for period detection
tomwalters@0: %
tomwalters@0: function [F0estSAI, MeanSAI, ModMeanSAI, F0estSAIprof, ProfSAI ] ...
tomwalters@0:              = CalF0estSAI(SAIval, fs, Frs, StrobeLoc, MaxF0, PreF0);
tomwalters@0: 
tomwalters@0: if nargin < 4, StrobeLoc = 5*fs/1000+1; end; % at 5 ms
tomwalters@0: if nargin < 5, MaxF0 = 400; end;
tomwalters@0: if nargin < 6, PreF0 = 150; end;
tomwalters@0: 
tomwalters@0: MinF0 = 80;
tomwalters@0: MeanSAI = mean(SAIval);
tomwalters@0: [NCh LenSAI] = size(SAIval);
tomwalters@0: 
tomwalters@0: [dummy nc1] = min(abs(Frs-500));
tomwalters@0: % wgtChSAI = [ (1:nc1)'/nc1; ones(NCh-nc1,1)]; useless
tomwalters@0: wgtChSAI = ones(NCh,1);
tomwalters@0: 
tomwalters@0: nTilt = StrobeLoc + fix(fs/max(PreF0,MinF0)); % F0 > 80 Hz
tomwalters@0: nnt = [0:(LenSAI-nTilt)];
tomwalters@0: wgtTISAI = [zeros(1,StrobeLoc-1), ones(1,nTilt-StrobeLoc),  ...
tomwalters@0: 		(1 - nnt/max(nnt)) ];
tomwalters@0: ModMeanSAI = mean((wgtChSAI*wgtTISAI).*SAIval);
tomwalters@0: 
tomwalters@0: %[dummy nc0] = min(abs(Frs-70)); % Frs > 100
tomwalters@0: nc0 = 1;
tomwalters@0: wgtProSAI = [zeros(1,nc0), (1 - ((nc0+1):nc1)/nc1) zeros(1,NCh-nc1)]';
tomwalters@0: 
tomwalters@0: ProfSAI = wgtProSAI.*mean(SAIval')';
tomwalters@0: [dummy mp] = max(ProfSAI);
tomwalters@0: F0estSAIprof = Frs(mp);
tomwalters@0: 
tomwalters@0: ThreshUV = 0.95;
tomwalters@0: 
tomwalters@0:    F0estSAI = 0;
tomwalters@0:    if  max(MeanSAI) == 0, F0estSAI = 0; return; end;
tomwalters@0: 
tomwalters@0:    [PeakVal1 PeakLoc1 ] = max(MeanSAI);
tomwalters@0:    if PeakLoc1 ~= StrobeLoc;	
tomwalters@0:   	disp([ 'StrobeLoc is strangely detected at ' int2str(PeakLoc1)]);
tomwalters@0: 	F0estSAI = 0;
tomwalters@0:  	return;
tomwalters@0:    end;
tomwalters@0: 
tomwalters@0:    ModMeanSAI2 = ModMeanSAI;
tomwalters@0:    nz = StrobeLoc+(0: fix(fs/MaxF0)-1);
tomwalters@0:    ModMeanSAI2(nz) = zeros(size(nz)); % suppress 1st peak
tomwalters@0:    [PeakVal2 PeakLoc2] = max(ModMeanSAI2);
tomwalters@0:    F0estSAI = fs/(PeakLoc2-StrobeLoc);
tomwalters@0: 	
tomwalters@0:    [TroughVal dummy] = min(ModMeanSAI(StrobeLoc:PeakLoc2));
tomwalters@0:    if TroughVal/PeakVal2 >= ThreshUV;
tomwalters@0: 	F0estSAI = 0;
tomwalters@0:    end; 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: