ls.m

《现代信号处理》中关于利用最小二乘法估计ARMA模型的参数

%--------------最小二乘法程序----------------------%
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')