function [l,threshold]=estimate(rr,S)
% 循環(huán)周期數(shù)S，可改為 20,40,80
k=fix(length(rr)/S);                                              % 相關(guān)延時(shí)，可改                                                                                
y=[zeros(1,k), rr];                                               
N=fix(length(rr)/S);                                              % 循環(huán)周期長度:2048,1024,512,204,102,51,25,12
q=fix(N);                                                         % 采樣點(diǎn)數(shù)，可改為Ｎ或Ｎ的整數(shù)分之一

for n=1:q  
    for  s=1:S 
        %Cxx(n,s) = xcorr(rr(n+(s-1)*N),'unbiased');
        Cxx(n,s) =(rr(n+(s-1)*N))*(rr(n+(s-1)*N));
    end  
        Cxx1(n) = sum(Cxx(n,:))/S;             
        Cxx_ave = mean(Cxx1(1,:));
        Cxx2 = Cxx1-Cxx_ave;  
  end

for m=1:q
     for  s=1:S
            %Cyy(m,s) = xcorr(y(m+(s-1)*N),'unbiased');
            Cyy(m,s) =(y(m+(s-1)*N))*(y(m+(s-1)*N));                     
     end 
            Cyy1(m) = sum(Cyy(m,:))/S;  
end   
       COV = xcov(Cxx1,Cyy1);  
       lj=fix(sqrt(length(COV)));
       li = length(COV(1:lj^2)); 
       COV=reshape(COV(1:lj^2),lj,lj);
       Cxx3=reshape(Cxx2(1:lj),1,lj);    
       l = S * Cxx3*inv(COV)*Cxx3'
      
       threshold=chi2inv(0.7,q);                               % 原算法的固定閾值