ar2.m

来自「这是一个关于MATLAB的源代码」· M 代码 · 共 29 行

function [a, e] = AR2(x, p)
% AR自相关算法谱估计，输入采样数据x，预测阶数p，输出模型传输函数系数a
% 返回H(z)=e/(1+A(z))中的A，e是sigma^2
N = length(x);
A = zeros(p, N+p);
for n=1:p
    for m=1:N+p
        if m<=n || m-n>=N+1
            A(n, m) = 0;
        else 
            A(n, m) = x(m-n);  %A'
        end
    end
end
xf = zeros(N+p, 1);
xf(1:N) = x';
a = -1*pinv(A')*xf;
e = (xf'*xf+a'*A*xf)/(N+p);
%test
% clear;
% close all;
% a = [1, -1.3817, 1.5632, -0.8843, 0.4096];
% b = [1, 0.3544, 0.3508, 0.1736, 0.2401];
% x = xx(a, b, 200, 256);
% [coa, coe] = AR2(x, 4);
% coa = [1;coa];
% cob = sqrt(coe);
% [w,P] = spectrum(coa, cob, 128);
% plot(w, P);