www.pudn.com > cyclostationary_toolbox.rar > cyclic_cross_covariance_fast.m, change:1998-06-24,size:3090b

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