[求助]循环自相关函数

function R=cyclic_autocorrelation(x,alpha,max_tau)
R=cyclic_cross_correlation(x,x,alpha,max_tau);


function R=cyclic_cross_correlation(x,y,alpha,max_tau)
%
% CYCLIC_CROSS_CORRELATION
%              
%              calculates the cyclic cross correlation between
%              two signals x,y at frequency alpha
%            
%              R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)}
%              for k=0 ... 2*pi/alpha-1
%
% USAGE
%              R=cyclic_cross_correlation(x,y,alpha,max_tau)
%
%              calculate cross correlation up to max_tau time lags

% File: cyclic_cross_correlation.m


% Simple error checks
if nargin~=4
  error('Incorrect number of arguments for function cyclic_cross_correlation');
end
if alpha>2*pi
  error('Cyclic frequency must be less than 2 pi in function cyclic_cross_correlation');
end


T=ceil(2*pi/alpha)-1;
lx=length(x);
t=0:lx-1;
R=zeros(max_tau*2+1,T+1);


% Compute even time shift segments
for tau=-max_tau:2:max_tau
  for k=0:T
    R(tau+1+max_tau,k+1)=mean(x(1:lx-max_tau-tau).*y(max_tau+tau+1:lx) ...
    .*exp(-j*k*alpha*t(1+(max_tau+tau)/2:lx-(max_tau+tau)/2)));
  end
end

% Compute odd time shift segments
t=t+0.5;
for tau=-max_tau+1:2:max_tau
  for k=0:T
    R(tau+1+max_tau,k+1)=mean(x(1:lx-tau-max_tau).*y(max_tau+tau+1:lx) ...
    .*exp(-j*k*alpha*t(1+(max_tau+tau-1)/2:lx-(max_tau+tau+1)/2)));
  end
end
这里的alpha似乎是分辨率，我这样子调用的
clear;
t=(0:1024-1)/1024;
x=sin(200*pi*t+30*pi/180).*(1+1.5*cos(30*pi*t));
R=cyclic_autocorrelation(x,1/1024,10);mesh(abs(R));
仅供参考，呵呵