所用函数一：

function [A,Sgm2] = AR(x,P)
% x - 输入数据，列向量
% P - AR 阶数,标量
% A - AR 系数[a1 a2 ...aP]，行向量;

A = zeros(1,P); % 初始化参数矩阵，行向量
Rx = Rxx(x,P); % 无偏自相关估计，size = P+1

% p = 1 时候的Yule-Walker方程解,注意，MATLAB下标从1开始
A(1) = - Rx(1+1)/Rx(0+1);
Sgm2 = Rx(0+1)*(1 - A(1)^2);

% p = 1->P 时候的递推,注意，MATLAB下标从1开始
for p = 1:P-1
k = 1:p;
K = -(Rx((p+1)+1) + A(k)*Rx((p+1-k)+1))/Sgm2;
Sgm2 = Sgm2 * (1 - K*K);

A(k) = A(k) + K * A(p+1-k);
A(p+1) = K;
end

所用函数二：

function Rx = Rxx(x,P)
% x - 自变量
% P - 自相关数目
% Rxx[m] = sum(x[n]*x[n+m]),n = 0->inf,m = 0->P;

N = length(x);
sizex = size(x);
if sizex(1)==1
Rx = zeros(1,P+1); % x 为行向量那么Rx也为行向量
for m = 0:P
Rx(m+1) = x(m+1:N)*x(1:N-m)';
end
elseif sizex(2)==1
Rx = zeros(P+1,1); % x 为列向量那么Rx也为列向量
for m = 0:P
Rx(m+1) = x(m+1:N)'*x(1:N-m);
end
end

所用函数三：

function [H,f] = ar2fhz(A,fs,df)
% A - AR 的系数[a1 a2 ... ap]，行向量;
% fs - 指定的频率范围
% dim - 分辨率,角度；
% H - 功率谱密度函数，行向量
% df - 频率分辨率，标量

P = length(A);

% Hz分辨率->弧度分辨率
dw = df/fs*2*pi;

w = 0:dw:pi;
z = exp(-j*w);
H = zeros(size(w));
for i = 1:P
H = H + A(i)*z.^i;
end
f = w/2/pi*fs;
H = 1./(1 + H);
H = abs(H).^2;

测试

tic
fs = 600; % 采样率
t = 0:1/fs:1;
f1 = 100; % 频率100赫兹
f2 = 160; % 频率160赫兹
f3 = 230; % 频率230赫兹

x1 = sin(2*pi*f1*t);
x2 = 1.2*sin(2*pi*f2*t);
x3 = 1.5*sin(2*pi*f3*t);
x = x1 + x2 + x3;

df = 1; % 频域间隔

[a,~] = AR(x',40);
[H,f] = ar2fhz(a,fs,df);
plot(f,log(H),'r')
toc

结果

还行，还能用。。。