求助关于循环谱的问题

施虹

看了好几天的资料，还是稀里糊涂，哪位大大能够用浅显的语言帮我说说，循环谱的实现过程啊？
这里有个循环谱FAM解法的程序，高手们帮忙分析一下，偶看不懂！:@L
function S=fft_accumulation_method(x,y,N,a,g,L)
%
% FFT_ACCUMULATION_METHOD
%              estimate the spectral correlation density using the FFT
%              accumulation method
%使用FAM估算循环谱
%              Reference:
%              Roberts R. S., Brown, W. A. and Loomis, H. H.
%              "Computationally efficient algorithms for cyclic
%              spectral analysis" IEEE Signal Processing Magazine
%              8(2) pp38-49 April 1991
%              
%              Input parameters are:
%              x,y signals
%              N   length of time window used for estimating frequency
%                  segments (should be power of 2)
%              a   window used for smoothing segments
%              g   window for smoothing correlation
%              L   decimation factor
%              
% USAGE        S=fft_accumulation_method(x,y,N,a,g,L)
%
%              e.g s=fft_accumulation_method(s1,s1,64,'hamming','hamming',1)

if ~exist('L')
  L=1;
end

lx=length(x);
if (length(y)~=lx)
  error('Time series x and y must be same length')
end


n=0:floor((lx-N)/L);
ln=length(n);
a=feval(a,N)';
g=feval(g,ln)';
g=g/sum(g);
a=a/sum(a);

X=zeros(N,ln);
Y=zeros(N,ln);

Ts=1/N;

for i=1:ln
  n_r=n(i)*L+(1:N);
  X(:,i)=fftshift(fft(a.*x(n_r)))';
  Y(:,i)=fftshift(conj(fft(a.*y(n_r))))';   
end


lnd2=floor(ln/2);
lnd4=floor(ln/4)+1;
ln3d4=lnd4+lnd2-2;
S=zeros(2*N-1,lnd2*(N+1));

for alpha=-N/2+1:N/2-1
  for f=-N/2:N/2-1
    f1=f+alpha;
    f2=f-alpha;
    if ((abs(f1)<N/2)&(abs(f2)<N/2))
      f1=f1+N/2;
      f2=f2+N/2;
      fsh=fftshift(fft(g.*X(f1,:).*Y(f2,:)));
      S(2*f+N,(alpha+N/2)*lnd2+(1:lnd2-1))=fsh(1,lnd4:ln3d4);
    end
    f1=f+alpha;
    f2=f-alpha+1;
    if ((abs(f1)<N/2)&(abs(f2)<N/2))
      f1=f1+N/2;
      f2=f2+N/2;
      fsh=fftshift(fft(g.*X(f1,:).*Y(f2,:)));
      S(2*f+N+1,(alpha+N/2)*lnd2+(1:lnd2-1)-lnd4+1)=fsh(1,lnd4:ln3d4);
    end
  end  
end

先说声谢谢了！