cyclic_3rd_order_cumulant.m

来自「计算循环平稳信号的3阶累积量」· M 代码 · 共 62 行

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