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

```function C3=cyclic_3rd_order_cumulant_fast(x1,x2,x3,T,max_tau)
%
% CYCLIC_3RD_ORDER_CUMULANT_FAST
%
%              calculates the cyclic third order cumulant of
%              three signals x1,x2,x3 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.
%
%              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_fast(x1,x2,x3,alpha,max_tau)
%

% File: cyclic_3rd_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~=5
error('Incorrect number of arguments for function cyclic_third_order_cumulant_fast');
end
if T(1)<1
error('Synchronous period must be larger than 1 in function cyclic_third_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

% Calculate synchronous averages from signals
mx1=synchronous_average(x1,T);
mx2=synchronous_average(x2,T);
mx3=synchronous_average(x3,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;
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))];
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=[];
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];
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;
end

% Compute time varying third order cumulant
r=zeros(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
temp(:)=x1(t).*x2(t+tau1).*x3(t+tau2);
r(tau1+1,tau2+1,:)=mean(temp');
end
end

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

```