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