cyclic_cross_covariance.m

循环自相关函数工具箱

function R=cyclic_cross_covariance(x,y,alpha,max_tau)
%
% CYCLIC_CROSS_COVARIANCE
%              
%              calculates the cyclic cross covariance 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 ... 1/alpha
%              With signals x and y having their periodic average removed
%
% USAGE
%              R=cyclic_cross_covariance(x,y,alpha,max_tau)
%
%              calculate cross covariance up to max_tau time lags

% File: cyclic_cross_covariance.m
% Last Revised: 23/4/98
% 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');
end
if alpha>2*pi
  error('Cyclic frequency must be less than 2 pi in function cyclic_cross_covariance');
end




% Remove cyclic mean from signals
cmx=cyclic_mean(x,alpha);
cmy=cyclic_mean(y,alpha);
lx=length(x);
t=0:lx-1;
T=ceil(2*pi/alpha)-1;
for k=1:lx
  x(k)=x(k)-1/(2*pi)*sum(cmx.*exp(j*alpha*(0:T)*(k-1)));
  y(k)=y(k)-1/(2*pi)*sum(cmy.*exp(j*alpha*(0:T)*(k-1)));
end

R=zeros(2*max_tau,T+1);


% Compute even time shift segments
for tau=-max_tau:2:max_tau
  for k=0:T
    R(tau+1+max_tau,k+1)=mean(x(1:lx-tau-max_tau).*y(tau+1+max_tau:lx) ...
       .*exp(j*k*alpha*t(1+(tau+max_tau)/2:lx-(tau+max_tau)/2)));
  end
end

% Compute odd time shift segments
t=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(tau+1+max_tau:lx) ...
       .*exp(j*k*alpha*t(1+(tau+max_tau-1)/2:lx-(tau+max_tau+1)/2)));
  end
end