cyclostationary_toolbox/cyclic_3rd_order_cumulant.m 0100644 0031654 0031654 00000003034 06436565164 0023325 0 ustar 00acmc acmc 0000304 0001726 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 checks
if nargin~=5
error('Incorrect number of arguments for function cyclic_3rd_order_cumulant');
end
if alpha>2*pi
error('Cyclic frequency must be less than 2 pi in function cyclic_3rd_order_cumulant');
end
% Remove cyclic mean from signals
cmx1=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)));
end
C3=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
end
end
cyclostationary_toolbox/cyclic_3rd_order_cumulant_fast.m 0100644 0031654 0031654 00000006074 06436555370 0024347 0 ustar 00acmc acmc 0000304 0001726 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
cyclostationary_toolbox/cyclic_4th_order_cumulant.m 0100644 0031654 0031654 00000004173 06436565216 0023337 0 ustar 00acmc acmc 0000304 0001726 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
