cyclic_cross_covariance_fast.m

循环自相关函数工具箱

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 checks
if nargin~=4
  error('Incorrect number of arguments for function cyclic_cross_covariance_fast');
end
if 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 shapes
if 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 signals
mx=synchronous_average(x,T);
my=synchronous_average(y,T);


% Remove excess samples due to non-integer sampling
% and renove cyclic mean from signal
if 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 correlation
r=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 shift
R=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