katkost@3: function f0 = f0detectionyin(x,fs,ws,minf0,maxf0) 
katkost@3: % x: input signal, fs: sampling rate, ws: integration window length, minf0: minimum f0, maxf0: maximum f0 
katkost@3: % f0: fundamental frequency detected in Hz 
katkost@3: maxlag = ws-2;                            % maximum lag 
katkost@3: th = 0.1;                                 % set threshold 
katkost@3: d = zeros(maxlag,1);                      % init variable 
katkost@3: d2 = zeros(maxlag,1);                     % init variable 
katkost@3: % compute d(tau) 
katkost@3: x1 = x(1:ws);  
katkost@3: cumsumx = sum(x1.^2); 
katkost@3: cumsumxl = cumsumx; 
katkost@3: xy = xcorr(x(1:ws*2),x1); 
katkost@3: xy = xy(ws*2+1:ws*3-2); 
katkost@3: for lag=0:maxlag-1 
katkost@3:     d(1+lag) = cumsumx + cumsumxl - 2*xy(1+lag); 
katkost@3:     cumsumxl = cumsumxl - x(1+lag).^2 + x(1+lag+ws+1)^2; 
katkost@3: end 
katkost@3: cumsum = 0; 
katkost@3: % compute d'(tau) 
katkost@3: d2(1) = 1; 
katkost@3: for lag=1:maxlag-1 
katkost@3:     cumsum = cumsum + d(1+lag); 
katkost@3:     d2(1+lag) = d(1+lag)*lag./cumsum; 
katkost@3: end 
katkost@3: % limit the search to the target range 
katkost@3: minf0lag = 1+round(fs./minf0);            % compute lag corresponding to minf0 
katkost@3: maxf0lag = 1+round(fs./maxf0);            % compute lag corresponding to maxf0 
katkost@3: if (maxf0lag>1 && maxf0lag<maxlag) 
katkost@3:     d2(1:maxf0lag) = 100;                 % avoid lags shorter than maxf0lag 
katkost@3: end 
katkost@3: if (minf0lag>1 && minf0lag<maxlag) 
katkost@3:         d2(minf0lag:end) = 100;               % avoid lags larger than minf0lag 
katkost@3: end 
katkost@3: % find the best candidate 
katkost@3: mloc = 1 + find((d2(2:end-1)<d2(3:end)).*(d2(2:end-1)<d2(1:end-2)));  % minima 
katkost@3: candf0lag = 0; 
katkost@3: if (length(mloc)>0) 
katkost@3:     I = find(d2(mloc)<th); 
katkost@3:     if (length(I)>0) 
katkost@3:         candf0lag = mloc(I(1)); 
katkost@3:     else 
katkost@3:         [Y,I2] = min(d2(mloc)); 
katkost@3:         candf0lag = mloc(I2); 
katkost@3:     end 
katkost@3:     candf0lag = candf0lag;                % this is zero-based indexing     
katkost@3:     if (candf0lag>1 & candf0lag<maxlag) 
katkost@3:         % parabolic interpolation 
katkost@3:         lval = d2(candf0lag-1);  
katkost@3:         val = d2(candf0lag);  
katkost@3:         rval= d2(candf0lag+1);  
katkost@3:         candf0lag = candf0lag + .5*(lval-rval)./(lval-2*val+rval);   
katkost@3:     end 
katkost@3: end         
katkost@3: ac = min(d2);  
katkost@3: f0lag = candf0lag-1;                      % convert to zero-based indexing 
katkost@3: f0 = fs./f0lag;                           % compute candidate frequency in Hz 
katkost@3: if (ac > 0.2)                             % voiced/unvoiced threshold 
katkost@3:     f0 = 0;                               % set to unvoiced 
katkost@3: end 
katkost@3: