ar.m

来自「这是关于功率谱的Matlab仿真程序」· M 代码 · 共 27 行

function [a, e] = AR(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
        A(n, m) = x(p-n+m);  %A'
    end
end
xf = x(p+1:N);
[V1, D1] = eig(A*A');
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] = AR(x, 4);
% coa = [1;coa];
% cob = sqrt(coe);
% [w,P] = spectrum(coa, cob, 128);
% plot(w, P);