www.pudn.com > cyclostationary_toolbox.rar > cyclic_cross_correlation_fast.m, change:1998-06-08,size:2767b

function R=cyclic_cross_correlation_fast(x,y,T,max_tau)
%
% CYCLIC_CROSS_CORRELATION_FAST
%
%              calculates the cyclic cross correlation 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
%
%              R=cyclic_cross_correlation_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_correlation_fast.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_correlation_fast');
end
if T(1)<1
error('Synchronous period must be larger than 1 in function cyclic_cross_correlation_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

% Remove excess samples due to non-integer sampling
if length(T)==1
% remove jitter samples if non-integer T
if T~=floor(T)
cp=1;np=1;
while cp+T<length(x)
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
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)];
ny=[ny  y(cp:cp+T-1)];
end
nx=[nx x(rot_positions(n+1):rot_positions(n+1)+2*max_tau+1)];
x=nx;
ny=[ny y(rot_positions(n+1):rot_positions(n+1)+2*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