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 Tif 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)];endtemp=zeros(floor(T),n);z=(1:floor(T)*n);% calculate cyclic mean - time domain averagetemp(:)=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;