function R=cyclic_cross_covariance_fast(x,y,T,max_tau)%% CYCLIC_CROSS_COVARIANCE_FAST%              %              calculates the cyclic cross covariance between%              two signals x,y at frequency alpha=1/T using a fast%              approximation based on the synchronous average of the%              time varying autocorrelation.  Fundamental signal%              period can be defined as a single period or as a sequence%              of once per period pulse times.%             %              R(k*alpha,tau)=E{ x(t-tau/2) y(t+tau/2)%                                exp(-jk(alpha)t)}%              for k=0 ... 1/alpha%              With signals x and y having their periodic average removed%% USAGE%              R=cyclic_cross_covariance_fast(x,y,T,max_tau)%%              calculate cross correlation up to max_tau time lags%%              if T is a scalar, then T defined the fundamental%              cyclic period%%              if T is a vector, then T defines a series of once%              per revolution impulses% File: cyclic_cross_covariance_fast.m% Last Revised: 23/4/97% 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_covariance_fast');endif T(1)<1  error('Synchronous period must be larger than 1 in function cyclic_cross_covariance_fast');end% Ensure that vectors are the right sizes and shapesif length(x)>length(y)  x=x(1:length(y));end[rows,cols]=size(x);if rows>cols  x=x';end[rows,cols]=size(y);if rows>cols  y=y';end% Calculate synchronous averages from signalsmx=synchronous_average(x,T);my=synchronous_average(y,T);% Remove excess samples due to non-integer sampling% and renove cyclic mean from signalif length(T)==1    cp=1;np=1;  while cp+T<length(x)    x(cp:cp+floor(T)-1)=x(cp:cp+floor(T)-1)-mx;    y(cp:cp+floor(T)-1)=y(cp:cp+floor(T)-1)-my;    cp=cp+floor(T);    np=np+T;    if (np-cp)>1      x=[x(1:cp-1) x(cp+1:length(x))];      y=[y(1:cp-1) y(cp+1:length(y))];      np=np-1;    end  end  n=floor((length(x)-2*max_tau-1)/T);else  % extract time series correlated with periodic pulses  rot_positions=T;  T=floor(median(diff(rot_positions)));  nx=[];  ny=[];  n=length(rot_positions)-2;  for k=1:n;    cp=rot_positions(k);    nx=[nx; x(cp:cp+T-1)-mx];    ny=[ny; y(cp:cp+T-1)-my];  end  nx=[nx;x(rot_positions(n+1):rot_positions(n+1)+max_tau+1)-mx(1:max_tau+1)];  x=nx;  ny=[ny;y(rot_positions(n+1):rot_positions(n+1)+max_tau+1)-my(1:max_tau+1)];  y=ny;end% Compute time varying cross correlationr=zeros(2*max_tau+1,floor(T));temp=zeros(floor(T),n);t=(1:floor(T)*n); for k=-max_tau:max_tau      temp(:)=x(t+max_tau).*y(t+k+max_tau);   r(k+1+max_tau,:)=mean(temp');end% Take DFT of time varying correlation with appropriate phase change% to compensate for time shiftR=zeros(2*max_tau+1,floor(T));for k=-max_tau:max_tau    R(k+1+max_tau,:)=exp(-i*pi*((0:floor(T)-1)/T)*k).*fft(r(k+1+max_tau,:))/T;end