cyclic_4th_order_cumulant_fast.m

循环自相关函数工具箱

function C4=cyclic_4th_order_cumulant_fast(x1,x2,x3,x4,T,max_tau)
%
% CYCLIC_4TH_ORDER_CUMULANT_FAST
%              
%              calculates the cyclic fourth order cumulant of 
%              four signals x1,x2,x3,x4 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.
%             
%              C4(k*alpha,tau1,tau2)=E{(x1(t)-E{x1(t)}) *
%                                      (x2(t+tau1)-E{x2(t+tau1)} *
%             			       (x3(t+tau2)-E{x3(t+tau2)} *
%             			       (x4(t+tau2)-E{x4(t+tau2)} *
%                                      exp(-jk(alpha)t)  }
%              for k=0 ... 1/alpha
%             
%
% USAGE
%              C4=cyclic_4th_order_cumulant_fast(x1,x2,x3,x4,alpha,max_tau)
%

% File: cyclic_4th_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~=6
  error('Incorrect number of arguments for function cyclic_4th_order_cumulant_fast');
end
if T(1)<1
  error('Synchronous period must be larger than 1 in function cyclic_4th_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
[rows,cols]=size(x4);
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);
mx4=synchronous_average(x4,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;
    x4(cp:cp+floor(T)-1)=x4(cp:cp+floor(T)-1)-mx4;
    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))];
      x4=[x4(1:cp-1);x4(cp+1:length(x4))];
      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=[];
  nx4=[];
  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];
    nx4=[nx4; x4(cp:cp+T-1)-mx4];
  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;  
  nx4=[nx4;x4(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx4(1:tau+1)];
  x4=nx4;
end


% Compute time varying third order cumulant
r=zeros(max_tau+1,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
    for tau3=0:max_tau
      temp(:)=x1(t).*x2(t+tau1).*x3(t+tau2).*x4(t+tau3);M4=mean(temp');
      temp(:)=x1(t).*x2(t+tau1);c2_12=mean(temp');
      temp(:)=x3(t+tau2).*x4(t+tau3);c2_34=mean(temp');
      temp(:)=x1(t).*x3(t+tau2);c2_13=mean(temp');
      temp(:)=x2(t+tau1).*x4(t+tau3);c2_24=mean(temp');
      temp(:)=x1(t).*x4(t+tau3);c2_14=mean(temp');
      temp(:)=x2(t+tau1).*x3(t+tau2);c2_23=mean(temp');     
      r(tau1+1,tau2+1,tau3+1,:)=M4-c2_12.*c2_34-c2_13.*c2_24-c2_14.*c2_23;
    end
  end
end


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