function fangzhen1_tls()

%生成增廣矩陣Re
M=95;
pe=30;
loops=20;
fvMatrix=zeros(2,loops);
arMatrix=zeros(4,loops);
kk=0;
for loop=1:1:loops
    while kk~=4
         Re=zeros(M,pe+1);
         %w=randn(2000,1);
         w=wgn(1,2000,0);
         for n=1:1:128
             x(n)=(20^(1/2))*sin(2*pi*0.2*n)+(2^(1/2))*sin(2*pi*0.213*n)+w(n+loop*50);
             %x(n)=(20^(1/2))*sin(2*pi*0.2*n);
         end
         Rxx=xcorr(x,'unbiased');
         for i=1:1:M
               for j=1:1:pe+1
               Re(i,j)=Rxx(pe+i+1-j+128);
               end
         end

         %求Re的有效秩kk
        [U,S,V]=svd(Re);
        bottom = 0;
        for hh=1:1:pe+1
            bottom = bottom + S(hh,hh)^2;
        end
        overhead=0;
        v = overhead/bottom;
        kk=0;
  
        while v<0.998
            kk = kk + 1;
            overhead = overhead + S(kk,kk)^2;
             v = overhead/bottom;
        end
    end     
    kk;
    %求Sp
    Sp=zeros(kk+1);
    for j=1:1:kk
         for i=1:1:pe+1-kk
             Vj=V(:,j);
             Sp=Sp+S(j,j)^2*(Vj([i:1:i+kk]))*(Vj([i:1:i+kk]))';
         end
    end
    invSp=inv(Sp);
    fc=invSp(:,1)/invSp(1,1);
    arMatrix(:,loop)=fc(2:5);
    fz=roots(fc);
    count=1;
    for j=1:2:kk
        fw(count)=atan(imag(fz(j))/real((fz(j))));
        fv(count)=abs(fw(count)/(2*pi));
        count = count+1;     
    end
    fv=sort(fv);
    fvMatrix(:,loop)=fv';
end
%分別求AR參數估計和頻率估計的均值和方差
arMatrix
ar_mean = mean(arMatrix')'
ar_var = var(arMatrix')'
fvMatrix
fv_mean = mean(fvMatrix')'
fv_var = var(fvMatrix')'