cyclic_cumulants_fast.m

循环自相关函数工具箱

function [m,v,s,k]=cyclic_cumulants_fast(x,T)
% 
% CYCLIC_CUMULANTS_FAST
%                  calculate synchronously averaged zero lag cumulants
%                  up to fourth order
%
% USGAE
%                  [m,v,s,k]=cyclic_cumulants_fast(x,T)
%
%                  where x is a cyclostationary time series
%                  and  T is the fundamental period of interest
%                  or a vector of period start positions
%               
%                  output parameters are:
%                  m - synchronous average ( cyclic mean )
%                  v - cyclic variance
%                  s - cyclic skewness
%                  k - cyclic kurosis

% File: cyclic_cumulants_fast.m
% Last Revised: 25/11/97
% Created: 25/11/97
% Author: Andrew C. McCormick
% (C) University of Strathclyde



%remove jitter samples for non-integer T

if length(T)==1
  cp=1;np=1;
  while cp+T<length(x)
     cp=cp+floor(T);
     np=np+T;
     if (np-cp)>1
        x=[x(1:cp-1);x(cp+1:length(x))];
        np=np-1;
     end
  end
  n=floor((length(x)-1)/T);
  nts=x;
else
  rot_positions=T;
  T=floor(median(diff(rot_positions)));
  nts=[];
  n=length(rot_positions)-2;
  for k=1:n;
    cp=rot_positions(k);
    nts=[nts; x(cp:cp+T-1)];
  end
  nts=[nts;x(rot_positions(n+1):rot_positions(n+1)+1)];
end

temp=zeros(floor(T),n);
z=(1:floor(T)*n);

% calculate cyclic mean - time domain average
temp(:)=nts(z);
m1=mean(temp');
m2=mean(temp'.^2);
m3=mean(temp'.^3);
m4=mean(temp'.^4);

m=m1;
v=m2-m1.^2;
s=(m3-3*m2.*m1+2*m1.^3)./m2.^(1.5);
k=(m4-3*m2.^2+4*m3.*m1+12*m2.*m1.^2-6*m1.^4)./m2.^2;