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

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

```