function C3=cyclic_3rd_order_cumulant(x1,x2,x3,alpha,max_tau)%% CYCLIC_3RD_ORDER_CUMULANT%              %              calculates the cyclic third order cumulant of %              three signals x1,x2,x3 at frequency alpha%             %              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(x,y,alpha,max_tau)%% File: cyclic_3rd_order_cumulant.m% Last Revised: 25/11/97% Created: 25/11/97% Author: Andrew C. McCormick% (C) University of Strathclyde% Simple error checksif nargin~=5  error('Incorrect number of arguments for function cyclic_3rd_order_cumulant');endif alpha>2*pi  error('Cyclic frequency must be less than 2 pi in function cyclic_3rd_order_cumulant');end% Remove cyclic mean from signalscmx1=cyclic_mean(x1,alpha);cmx2=cyclic_mean(x2,alpha);cmx3=cyclic_mean(x3,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)));endC3=zeros(max_tau,max_tau,T+1);ix=1:lx-max_tau-1;for tau1=0:max_tau  for tau2=0:max_tau    for k=0:T      C3(tau1+1,tau2+1,k+1)=mean(x1(ix).*x2(tau1+ix) ...	  .*x3(tau2+ix).*exp(j*k*alpha*t(ix)));    end  endend