www.pudn.com > cyclostationary_toolbox.rar > frac_delay_cyclic_ac.m, change:1998-05-05,size:2094b

```function R=frac_delay_cyclic_ac(x,N,alpha,max_tau)
%
% FRAC_DELAY_CYCLIC_AC
%
%              calculates the cyclic auto correlation for signal x
%              at frequency alpha resampling 1/alpha to N samples
%              max_tau is expressed at original sampling fequency
%
%              R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)}
%              for k=0 ... N-1
%
% USAGE
%              R=frac_delay_cyclic_ac(x,N,alpha,max_tau)
%
%              calculate cyclic autocorrelation up to max_tau time lags

% File: frac_delay_cyclic_ac.m
% Last Revised: 5/5/98
% Created: 5/5/98
% Author: Andrew C. McCormick
% (C) University of Strathclyde

lx=length(x);

Rxx=zeros(lx-2*max_tau,2*max_tau+1);

for k=-max_tau:max_tau
Rxx(:,k+max_tau+1)=[x(max_tau+1:lx-max_tau).*x(k+(max_tau+1:lx-max_tau))]';
end

r=zeros(N,2*max_tau+1);
w=zeros(N,1);

T=1/alpha;
for k=1:lx-2*max_tau
cycle=floor(k/T);
position=(k-cycle*T)/T*N+1;

fp=floor(position);cp=ceil(position);
if (fp==cp)
if (fp==0)
w(N)=w(N)+1;
r(N,:)=r(N,:)+Rxx(k,:);
else
w(fp)=w(fp)+1;
r(fp,:)=r(fp,:)+Rxx(k,:);
end
else
if (fp==0)
w(N)=w(N)+position-fp;
r(N,:)=r(N,:)+Rxx(k,:)*(position-fp);
else
w(fp)=w(fp)+position-fp;
r(fp,:)=r(fp,:)+Rxx(k,:)*(position-fp);
end
if cp>N
w(1)=w(1)+cp-position;
r(1,:)=r(1,:)+Rxx(k,:)*(cp-position);
else
w(cp)=w(cp)+cp-position;
r(cp,:)=r(cp,:)+Rxx(k,:)*(cp-position);
end
end
end

% If insufficient samples to estimate r, fill in blanks with
% the neighbouring sample in r.
if (sum(w==0)>0)
warning('Insufficient samples for desired resolution');
for k=1:N
if (w(k)==0)
if k>1
r(k,:)=r(k-1,:);
end
else
r(k,:)=r(k,:)/w(k);
end
end
else
r=r./(w*ones(1,2*max_tau+1));
end

% Take DFT of time varying correlation with appropriate phase change
% to compensate for time shift
R=zeros(2*max_tau+1,N);
for k=-max_tau:max_tau
R(k+1+max_tau,:)=exp(-i*pi*((0:N-1)/N)*k).*fft(r(:,k+1+max_tau)');
end

```