cyclic_3rd_order_cumulant_fast.m

循环自相关函数工具箱

function C3=cyclic_3rd_order_cumulant_fast(x1,x2,x3,T,max_tau)
%
% CYCLIC_3RD_ORDER_CUMULANT_FAST
%              
%              calculates the cyclic third order cumulant of 
%              three signals x1,x2,x3 at frequency alpha 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.
%             
%              C3(k*alpha,tau1,tau2)=E{(x1(t)-E{x1(t)}) *
%                                      (x2(t+tau1)-E{x2(t+tau1)} *
%             			       (x3(t+tau2)-E{x3(t+tau2)} *
%                                      exp(-jk(alpha)t)  }
%              for k=0 ... 1/alpha
%             
%
% USAGE
%              C3=cyclic_3rd_order_cumulant_fast(x1,x2,x3,alpha,max_tau)
%

% File: cyclic_3rd_order_cumulant_fast.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde


% Simple error checks
if nargin~=5
  error('Incorrect number of arguments for function cyclic_third_order_cumulant_fast');
end
if T(1)<1
  error('Synchronous period must be larger than 1 in function cyclic_third_order_cumulant_fast');
end

[rows,cols]=size(x1);
if rows>cols
  x1=x1';
end
[rows,cols]=size(x2);
if rows>cols
  x2=x2';
end
[rows,cols]=size(x3);
if rows>cols
  x3=x3';
end

% Calculate synchronous averages from signals
mx1=synchronous_average(x1,T);
mx2=synchronous_average(x2,T);
mx3=synchronous_average(x3,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(x1)
    x1(cp:cp+floor(T)-1)=x1(cp:cp+floor(T)-1)-mx1;
    x2(cp:cp+floor(T)-1)=x2(cp:cp+floor(T)-1)-mx2;
    x3(cp:cp+floor(T)-1)=x3(cp:cp+floor(T)-1)-mx3;
    cp=cp+floor(T);
    np=np+T;
    if (np-cp)>1
      x1=[x1(1:cp-1);x1(cp+1:length(x1))];
      x2=[x2(1:cp-1);x2(cp+1:length(x2))];
      x3=[x3(1:cp-1);x3(cp+1:length(x3))];
      np=np-1;
    end
  end
  n=floor((length(x1)-2*max_tau-1)/T);
else
  % extract time series correlated with periodic pulses
  rot_positions=T;
  T=floor(median(diff(rot_positions)));
  nx1=[];
  nx2=[];
  nx3=[];
  n=length(rot_positions)-2;
  for k=1:n;
    cp=rot_positions(k);
    nx1=[nx1; x1(cp:cp+T-1)-mx1];
    nx2=[nx2; x2(cp:cp+T-1)-mx2];
    nx3=[nx3; x3(cp:cp+T-1)-mx3];
  end
  nx1=[nx1;x1(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx1(1:tau+1)];
  x1=nx1;
  nx2=[nx2;x2(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx2(1:tau+1)];
  x2=nx2;
  nx3=[nx3;x3(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx3(1:tau+1)];
  x3=nx3;
end


% Compute time varying third order cumulant
r=zeros(max_tau+1,max_tau+1,floor(T));
temp=zeros(floor(T),n);
t=(1:floor(T)*n); 
for tau1=0:max_tau
  for tau2=0:max_tau
    temp(:)=x1(t).*x2(t+tau1).*x3(t+tau2);
    r(tau1+1,tau2+1,:)=mean(temp');
  end
end


% Take DFT of time varying toc
C3=zeros(max_tau+1,max_tau+1,floor(T));
for tau1=0:max_tau
  for tau2=0:max_tau
    C3(tau1+1,tau2+1,:)=fft(r(tau1+1,tau2+1,:))/T;
  end
end