cyclic_cross_correlation_fast.m

一些循环平稳子程序

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