frac_delay_cyclic_ac.m

循环自相关函数工具箱

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 Strathclyde


lx=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))]';
end

r=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
  end
end


% 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
  end
else
  r=r./(w*ones(1,2*max_tau+1));
end


% Take DFT of time varying correlation with appropriate phase change
% to compensate for time shift
R=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