function R=frac_delay_cyclic_ac(x,N,alpha,max_tau)%% FRAC_DELAY_CYCLIC_AC%              %              calculates the cyclic auto correlation for signal x%              at frequency alpha resampling 1/alpha to N samples%              max_tau is expressed at original sampling fequency%             %              R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)}%              for k=0 ... N-1%% USAGE%              R=frac_delay_cyclic_ac(x,N,alpha,max_tau)%%              calculate cyclic autocorrelation up to max_tau time lags% File: frac_delay_cyclic_ac.m% Last Revised: 5/5/98% Created: 5/5/98% Author: Andrew C. McCormick% (C) University of Strathclydelx=length(x);Rxx=zeros(lx-2*max_tau,2*max_tau+1);for k=-max_tau:max_tau  Rxx(:,k+max_tau+1)=[x(max_tau+1:lx-max_tau).*x(k+(max_tau+1:lx-max_tau))]';endr=zeros(N,2*max_tau+1);w=zeros(N,1);T=1/alpha;for k=1:lx-2*max_tau  cycle=floor(k/T);  position=(k-cycle*T)/T*N+1;  fp=floor(position);cp=ceil(position);  if (fp==cp)    if (fp==0)      w(N)=w(N)+1;      r(N,:)=r(N,:)+Rxx(k,:);    else      w(fp)=w(fp)+1;      r(fp,:)=r(fp,:)+Rxx(k,:);    end  else    if (fp==0)      w(N)=w(N)+position-fp;      r(N,:)=r(N,:)+Rxx(k,:)*(position-fp);    else      w(fp)=w(fp)+position-fp;      r(fp,:)=r(fp,:)+Rxx(k,:)*(position-fp);    end    if cp>N      w(1)=w(1)+cp-position;      r(1,:)=r(1,:)+Rxx(k,:)*(cp-position);    else      w(cp)=w(cp)+cp-position;      r(cp,:)=r(cp,:)+Rxx(k,:)*(cp-position);    end  endend% If insufficient samples to estimate r, fill in blanks with% the neighbouring sample in r.if (sum(w==0)>0)  warning('Insufficient samples for desired resolution');  for k=1:N    if (w(k)==0)      if k>1	r(k,:)=r(k-1,:);      end    else      r(k,:)=r(k,:)/w(k);    end  endelse  r=r./(w*ones(1,2*max_tau+1));end% Take DFT of time varying correlation with appropriate phase change% to compensate for time shiftR=zeros(2*max_tau+1,N);for k=-max_tau:max_tau    R(k+1+max_tau,:)=exp(-i*pi*((0:N-1)/N)*k).*fft(r(:,k+1+max_tau)');end