clc
total_len=100000;
L=16;
M=584;
Coff=48;
Po=0.01;
Fk=[-0.0253,-0.0342,-0.0358,-0.0167,0.0216,0.0649,0.0910,0.0819,0.0371,-0.0220,-0.0607,-0.0512,0.0079,0.0844,0.1269,0.0945,-0.0128,-0.1435,-0.2118,-0.1405,0.0946,0.4414,0.7859,1.0,1.0,0.7859,0.4414,0.0946,-0.1405,-0.2118,-0.1435,-0.0128,0.0945,0.1269,0.0844,0.0079,-0.0512,-0.0607,-0.0220,0.0371,0.0819,0.0910,0.0649,0.0216,-0.0167,-0.0358,-0.0342,-0.0253];
 %求Hemiian matrix Q,H;
  H=zeros(L,L+Coff-1);
  for i5=1:L
      for i6=1:Coff
        H(i5,i6+i5-1)=Fk(i6);
      end
   end
   G=zeros(L,L);
   Q=zeros(L,L);
   Qinv=zeros(L,L);
   G=H*H';
   Q=sqrt(G);
   Qinv=inv(Q);
   
%計算門限比 r1
 r1=(1+(L-1)*sqrt(2/(M*pi)))/(1+sqrt(2)*erfinv(1-2*(1-Po))*sqrt(2/M))
 r2=(L+M+1)/(M+2+sqrt(2)*erfinv(1-2*(1-Po))*sqrt(8*M+40+48/M))
  %求矩陣X(n); 
  a=zeros(1,total_len);
for i=1:total_len
    a(1,i)=simout(i);
end
% xn=zeros(1,M);
% for i1=1:M
%     xn(1,i1)=a(1,i1);
% end
y=zeros(1,L);
for i2=1:L
    for i3=1:M
       y(1,i2)= y(1,i2)+a(1,i3)*a(1,i3+i2-1)';
   end
end
 Rx=zeros(L,L);
 for i4=1:L
     for i5=1:L
         if i4==i5
             Rx(i4,i5)=y(1,1);
         else if i4<i5
                  Rx(i4,i5)=y(1,i5-i4+1);
              else
                   Rx(i4,i5)=y(1,i4-i5+1)';
               end
           end
       end
   end
end
 Rx=Rx/M;
    
% a=zeros(1,total_len);
% for i=1:total_len
%     a(1,i)=simout(i);
% end

% xn=zeros(L,M);

%for j=0:999
%    for i1=1:L       
%        for i2=1:M
%            xn(i1,i2)=a(1,i1+i2-1);
%        end
%    end
   
%    %求矩陣X(n)的協方差矩陣，Rx；
%    xc=xn';
%    Rx1=xn*xc;
%   Rx=zeros(L,L);
%   for i3=1:L
%       for i4=1:L
%           Rx(i3,i4)=norm(Rx1(i3,i4));
%       end
%   end

 Rxn=Qinv*Rx*Qinv;
 Tn1=0;
 Tn2=0;
 for i6=1:L
     for i8=1:L
         Tn1=Tn1+(norm(Rxn(i6,i8)))^2;
     end
 end
 for i7=1:L
     for i8=1:L
         if abs(i7-i8)<1
       Tn2=Tn2+(norm(Rxn(i7,i8)))^2;
         end
      end
 end
k=Tn1/Tn2