www.pudn.com > cyclostationary_toolbox.rar > cyclic_4th_order_cumulant.m, change:1997-11-25,size:2171b

```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

```