%--------------最小二乘法程序----------------------%
clear
% loops = 20;        % loops 是獨立實驗次數
% for loop=1:loops
   % 生成觀測數據
   n = [1:128];
   x = sqrt(20)*sin(2*pi*0.2*n) + sqrt(2)*sin(2*pi*0.213*n)+randn(1,128); % x(1,128)的行向量
   %w=randn(2000,1);
   %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);   %產生觀測數據的另一種方法，隨機性更好
   %end
   Rxx = xcorr(x,'unbiased');    % 相關矩陣
   % AR階數p給定為4或6
   p = 4;
   % p = 6
   q = p;                      % 用于諧波恢復的ARMA模型是特殊的APRM
   for i = 1:p
       for j = 1:p+1
           Re(i,j) = Rxx(p+i+1-j);
       end
   end                         % Re是修正方程的增廣矩陣
   %最小二乘法估計ARMA模型的參數
   A = Re(:,[2:p+1]);
   b = (-1)*Re(:,1);           % 修正方程 A*a = b     注意：a=[a1,a2,...ap]
   a = inv(A'*A)*A'*b
   % ARMA_LS(:,loop) = a;
   % 構造特征多項式A(z)
   Az(1) = 1;
   for i = 2:p+1
       Az(i) = a(i-1);
   end                         
   z = roots(Az)
   num = 1;
   for i = 1:2:p
       f_ls(num) = atan(imag(z(i))/real(z(i)));
       f_ls(num) = abs(f_ls(num))/(2*pi);
       num = num +1;
   end
   f_ls = sort(f_ls)
   %　Fz(:,loop) = f_ls;
   %end
% ARMA_LS
% ARMA_LS_mean = mean(ARMA_LS,2)
% ARMA_LS_std  = std(ARMA_LS')
% Fz
% Fz_mean = mean(Fz,2)
% Fz_std  = std(Fz')