cyclic_4th_order_cumulant.m

循环自相关函数工具箱

function C4=cyclic_4th_order_cumulant(x1,x2,x3,x4,alpha,max_tau)
%
% CYCLIC_4TH_ORDER_CUMULANT
%              
%              calculates the cyclic fourth order cumulant of 
%              three signals x1,x2,x3,x4 at frequency alpha
%             
%              C3(k*alpha,tau1,tau2,tau3)=E{(x1(t)-E{x1(t)}) *
%                                      (x2(t+tau1)-E{x2(t+tau1)} *
%             			       (x3(t+tau2)-E{x3(t+tau2)} *
%             			       (x4(t+tau3)-E{x4(t+tau3)} *
%                                      exp(-jk(alpha)t)  }
%              for k=0 ... 1/alpha
%             
%
% USAGE
%              C3=cyclic_4th_order_cumulant(x,y,alpha,max_tau)
%

% File: cyclic_4th_order_cumulant.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');
end
if alpha>2*pi
  error('Cyclic frequency must be less than 2 pi in function cyclic_4th_order_cumulant');
end


% Remove cyclic mean from signals
cmx1=cyclic_mean(x1,alpha);
cmx2=cyclic_mean(x2,alpha);
cmx3=cyclic_mean(x3,alpha);
cmx4=cyclic_mean(x4,alpha);
lx=length(x1);
t=0:lx-1;
T=ceil(2*pi/alpha)-1;
for k=1:lx
  x1(k)=x1(k)-1/(2*pi)*sum(cmx1.*exp(j*alpha*(0:T)*(k-1)));
  x2(k)=x2(k)-1/(2*pi)*sum(cmx2.*exp(j*alpha*(0:T)*(k-1)));
  x3(k)=x3(k)-1/(2*pi)*sum(cmx3.*exp(j*alpha*(0:T)*(k-1)));
  x4(k)=x4(k)-1/(2*pi)*sum(cmx4.*exp(j*alpha*(0:T)*(k-1)));
end

C4=zeros(max_tau,max_tau,max_tau,T+1);

ix=1:lx-max_tau-1;

for tau1=0:max_tau
  for tau2=0:max_tau
    for tau3=0:max_tau
      for k=0:T
	M4=mean(x1(ix).*x2(tau1+ix).*x3(tau2+ix) ...
	    .*x4(tau3+ix).*exp(j*k*alpha*t(ix)));
	c2_12=mean(x1(ix).*x2(tau1+ix).*exp(j*k*alpha*t(ix)));
	c2_34=mean(x3(tau2+ix).*x4(tau3+ix).*exp(j*k*alpha*t(ix)));
	c2_13=mean(x1(ix).*x3(tau2+ix).*exp(j*k*alpha*t(ix)));
	c2_24=mean(x2(tau1+ix).*x4(tau3+ix).*exp(j*k*alpha*t(ix)));
	c2_14=mean(x1(ix).*x4(tau3+ix).*exp(j*k*alpha*t(ix)));
	c2_23=mean(x2(tau1+ix).*x3(tau2+ix).*exp(j*k*alpha*t(ix)));
	C4(tau1+1,tau2+1,tau3+1,k+1)=M4-c2_12*c2_34-c2_13*c2_24-c2_14*c2_23;
      end
    end
  end
end