svd-tls_arma.m

来自「一种用总体最小二乘法来估计ARMA模型的实例」· M 代码 · 共 66 行

%ARMA谱估计
clear all
%设置观测数据
N1=128;
n=(0:N1-1);
x1=sqrt(20).*sin(2.*pi.*0.2*n)+sqrt(2).*sin(2.*pi.*0.213*n);
x=x1+randn(1,N1);
%设置数据长度
qe=100;
M=80;
pe=50;
%计算信号的自相关函数
rx=xcorr(x,'unbiased');
%构造样本相关函数矩阵Re
for j=1:pe+1,
     for i=1:M,
        
          Re(i,j)=rx(qe+i-j+1);
   end
end
q=qe;
%误差向量b
m=1:M;
for i=q+1:q+M,
    b(m)=rx(i);
end
%构造增广矩阵
B=[b',Re];
%用SVD法计算阶数p
[U,D,V]=svd(B);
for i=1:min(M,pe+1),
    j=i;
    dd=D(i,j)/D(1,1);
    if dd<0.05,
        break
    end
end
p=i-1;
R=rank(B);
%用TLS算法估计参数a
sp=0;
for j=1:p,
    for i=1:pe+1+1-p,
        vij=V(i:i+p,j);
        d(j)=D(j,j);
        s=d(j)^2*vij*vij';
        sp=sp+s;
    end
end
sp1=inv(sp);
for i=1:p,
    a(i)=sp1(i+1,1)/sp1(1,1);
end
A=[1,a];
figure;
freqz(1,A,1024,1);