Dawn@4: function r=yink(p,fileinfo) 
Dawn@4: % YINK - fundamental frequency estimator
Dawn@4: % new version (feb 2003)
Dawn@4: %
Dawn@4: %
Dawn@4: 
Dawn@4: %global jj;
Dawn@4: %jj=0;
Dawn@4: % process signal a chunk at a time
Dawn@4: idx=p.range(1)-1;
Dawn@4: totalhops=round((p.range(2)-p.range(1)+1) / p.hop);
Dawn@4: r1=nan*zeros(1,totalhops);r2=nan*zeros(1,totalhops);
Dawn@4: r3=nan*zeros(1,totalhops);r4=nan*zeros(1,totalhops);
Dawn@4: idx2=0+round(p.wsize/2/p.hop);
Dawn@4: while (1)
Dawn@4: 	start = idx+1;
Dawn@4: 	stop = idx+p.bufsize;
Dawn@4: 	stop=min(stop, p.range(2));
Dawn@4: 	xx=sf_wave(fileinfo, [start, stop], []);
Dawn@4: % 	if size(xx,1) == 1; xx=xx'; end
Dawn@4: 	xx=xx(:,1); 			% first channel if multichannel
Dawn@4: 	[prd,ap0,ap,pwr]=yin_helper(xx,p);
Dawn@4: 	n=size(prd ,2);
Dawn@4: 	if (~n) break; end;
Dawn@4: 	idx=idx+n*p.hop;
Dawn@4: 
Dawn@4: 	r1(idx2+1:idx2+n)= prd;	
Dawn@4: 	r2(idx2+1:idx2+n)= ap0;
Dawn@4: 	r3(idx2+1:idx2+n)= ap;
Dawn@4: 	r4(idx2+1:idx2+n)= pwr;
Dawn@4: 	idx2=idx2+n;
Dawn@4: end
Dawn@4: r.r1=r1; % period estimate
Dawn@4: r.r2=r2; % gross aperiodicity measure
Dawn@4: r.r3=r3; % fine aperiodicity measure
Dawn@4: r.r4=r4; % power
Dawn@4: sf_cleanup(fileinfo);
Dawn@4: % end of program
Dawn@4: 
Dawn@4: 
Dawn@4: 
Dawn@4: % Estimate F0 of a chunk of signal
Dawn@4: function [prd,ap0,ap,pwr]=yin_helper(x,p,dd)
Dawn@4: smooth=ceil(p.sr/p.lpf);
Dawn@4: x=rsmooth(x,smooth); 	% light low-pass smoothing
Dawn@4: x=x(smooth:end-smooth+1);
Dawn@4: [m,n]=size(x);
Dawn@4: maxlag = ceil(p.sr/p.minf0); 	
Dawn@4: minlag = floor(p.sr/p.maxf0);
Dawn@4: mxlg = maxlag+2; % +2 to improve interp near maxlag
Dawn@4: 
Dawn@4: 
Dawn@4: hops=floor((m-mxlg-p.wsize)/p.hop);
Dawn@4: prd=zeros(1,hops);
Dawn@4: ap0=zeros(1,hops);
Dawn@4: ap=zeros(1,hops);
Dawn@4: pwr=zeros(1,hops);
Dawn@4: if hops<1; return; end
Dawn@4: 
Dawn@4: % difference function matrix
Dawn@4: dd=zeros(floor((m-mxlg-p.hop)/p.hop),mxlg);
Dawn@4: if p.shift == 0		 % windows shift both ways
Dawn@4: 	lags1=round(mxlg/2) + round((0:mxlg-1)/2); 
Dawn@4: 	lags2=round(mxlg/2) - round((1:mxlg)/2);
Dawn@4: 	lags=[lags1; lags2];
Dawn@4: elseif p.shift == 1 % one window fixed, other shifts right
Dawn@4: 	lags=[zeros(1,mxlg); 1:mxlg]; 
Dawn@4: elseif p.shift == -1 % one window fixed, other shifts right
Dawn@4: 	lags=[mxlg-1:-1:0; mxlg*ones(1,mxlg)]; 
Dawn@4: else
Dawn@4: 	error (['unexpected shift flag: ', num2str(p.shift)]);
Dawn@4: end
Dawn@4: rdiff_inplace(x,x,dd,lags,p.hop);
Dawn@4: rsum_inplace(dd,round(p.wsize/p.hop));
Dawn@4: dd=dd';
Dawn@4: [dd,ddx]=minparabolic(dd); 	% parabolic interpolation near min
Dawn@4: cumnorm_inplace(dd);;		% cumulative mean-normalize
Dawn@4: 
Dawn@4: % first period estimate
Dawn@4: %global jj;
Dawn@4: for j=1:hops
Dawn@4: 	d=dd(:,j);
Dawn@4: 	if p.relflag
Dawn@4: 		pd=dftoperiod2(d,[minlag,maxlag],p.thresh); 
Dawn@4: 	else
Dawn@4: 		pd=dftoperiod(d,[minlag,maxlag],p.thresh); 
Dawn@4: 	end
Dawn@4: 	ap0(j)=d(pd+1);
Dawn@4: 	prd(j)=pd;
Dawn@4: end
Dawn@4: 
Dawn@4: % replace each estimate by best estimate in range
Dawn@4: range = 2*round(maxlag/p.hop);
Dawn@4: if hops>1; prd=prd(mininrange(ap0,range*ones(1,hops))); end
Dawn@4: %prd=prd(mininrange(ap0,prd)); 	
Dawn@4: 
Dawn@4: 
Dawn@4: % refine estimate by constraining search to vicinity of best local estimate
Dawn@4: margin1=0.6; 		
Dawn@4: margin2=1.8; 		
Dawn@4: for j=1:hops
Dawn@4: 	d=dd(:,j);
Dawn@4: 	dx=ddx(:,j);
Dawn@4: 	pd=prd(j);
Dawn@4: 	lo=floor(pd*margin1); lo=max(minlag,lo);
Dawn@4: 	hi=ceil(pd*margin2); hi=min(maxlag,hi);
Dawn@4: 	pd=dftoperiod(d,[lo,hi],0); 
Dawn@4: 	ap0(j)=d(pd+1);
Dawn@4: 	pd=pd+dx(pd+1)+1;	% fine tune based on parabolic interpolation
Dawn@4: 	prd(j)=pd;
Dawn@4: 
Dawn@4: 	% power estimates
Dawn@4: 	idx=(j-1)*p.hop;
Dawn@4: 	k=(1:ceil(pd))';
Dawn@4: 	x1=x(k+idx);
Dawn@4: 	x2=k+idx+pd-1;
Dawn@4: 	interp_inplace(x,x2);
Dawn@4: 	x3=x2-x1;
Dawn@4: 	x4=x2+x1;
Dawn@4: 	x1=x1.^2; rsum_inplace(x1,pd);
Dawn@4: 	x3=x3.^2; rsum_inplace(x3,pd);
Dawn@4: 	x4=x4.^2; rsum_inplace(x4,pd);
Dawn@4: 	
Dawn@4: 	x1=x1(1)/pd;
Dawn@4: 	x2=x2(1)/pd;
Dawn@4: 	x3=x3(1)/pd;
Dawn@4: 	x4=x4(1)/pd;
Dawn@4: 	% total power
Dawn@4: 	pwr(j)=x1;
Dawn@4: 
Dawn@4: 	% fine aperiodicity
Dawn@4: 	ap(j)=(eps+x3)/(eps+(x3+x4)); % accurate, only for valid min
Dawn@4: 
Dawn@4: 	%ap(j)
Dawn@4: 	%plot(min(1, d)); pause
Dawn@4: 
Dawn@4: 	prd(j)=pd;
Dawn@4: end
Dawn@4: 
Dawn@4: 
Dawn@4: 
Dawn@4: %cumulative mean-normalize 
Dawn@4: function y=cumnorm(x)
Dawn@4: [m,n]=size(x);
Dawn@4: y = cumsum(x);
Dawn@4: y = (y)./ (eps+repmat((1:m)',1,n)); % cumulative mean
Dawn@4: y = (eps+x) ./  (eps+y);