function R=cyclic_cross_correlation(x,y,alpha,max_tau)%% CYCLIC_CROSS_CORRELATION%              %              calculates the cyclic cross correlation between%              two signals x,y at frequency alpha%             %              R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)}%              for k=0 ... 2*pi/alpha-1%% USAGE%              R=cyclic_cross_correlation(x,y,alpha,max_tau)%%              calculate cross correlation up to max_tau time lags% File: cyclic_cross_correlation.m% Last Revised: 23/4/98% Created: 24/11/97% Author: Andrew C. McCormick% (C) University of Strathclyde% Simple error checksif nargin~=4  error('Incorrect number of arguments for function cyclic_cross_correlation');endif alpha>2*pi  error('Cyclic frequency must be less than 2 pi in function cyclic_cross_correlation');endT=ceil(2*pi/alpha)-1;lx=length(x);t=0:lx-1;R=zeros(max_tau*2+1,T+1);% Compute even time shift segmentsfor tau=-max_tau:2:max_tau  for k=0:T    R(tau+1+max_tau,k+1)=mean(x(1:lx-max_tau-tau).*y(max_tau+tau+1:lx) ...	.*exp(-j*k*alpha*t(1+(max_tau+tau)/2:lx-(max_tau+tau)/2)));  endend% Compute odd time shift segmentst=t+0.5;for tau=-max_tau+1:2:max_tau  for k=0:T    R(tau+1+max_tau,k+1)=mean(x(1:lx-tau-max_tau).*y(max_tau+tau+1:lx) ...	.*exp(-j*k*alpha*t(1+(max_tau+tau-1)/2:lx-(max_tau+tau+1)/2)));  endend