现代功率谱估计(2):Levinson-Durbin递推方法求解AR模型参数
p阶AR模型的差分方程形式和系统函数分别为:


令,则AR模型输出的功率谱密度为:

AR模型的系统输出信号为,计算输出信号的自相关函数:

其中最后结果的后面一项,输出信号和输入系统的白噪声的互相关函数满足关系:

其中为输入白噪声信号的方差,即噪声功率
因此上面的输入信号的自相关函数的最终取值应该分为两种情况讨论,分别式当中表示时间差的取和时,将代入,即:

写成矩阵形式的Yuler-Walker方程为:

其中实信号的自相关函数满足偶对称,即, 则上面的右边的自相关矩阵为Toeplitz矩阵,在确定AR模型阶数为时,通过求解AR模型的个系数,以及此时AR模型的输入白噪声的功率则可确定输出信号的功率谱密度分布。
p阶AR模型的系数分别为,其中的下标处第一个表示此时AR模型的阶数,也就是当前阶数情况下总共有几个系数,第二个下标i表示当前系数的序列索引。
以p=1,1阶AR模型为例,则写成矩阵形式的Yuler-Walker方程为:

求解未知数,为:

递推到p=2,2阶AR模型情况,写成矩阵形式的Yuler-Walker方程为:

求解p=2,2阶情况下的AR模型系数:

2阶AR模型的系数均可通过之前1阶AR模型的系数和估计出的自相关函数得到。
以此类推,k阶AR模型系数,白噪声功率均可通过第k-1阶AR模型的相关系数得到:

通过不断的往下一阶递推,可以得到较高阶AR模型,流程如下:

在循环递推开始之前,需要提前估计自相关函数,以及1阶情况的AR模型的相关系数和白噪声功率。
代码如下:
% levinson-durbin递推求解Yuler_Walker方程
clc; clear; close all;
N = 256;
p = 100;
n = 1: N;
fs = 2000;
t = n * 1/ fs;
f1 = 100; f2 = 230;
xn = 2 * sin(2 * pi * f1 * t + pi / 3) + sin(2 * pi * f2 * t + pi / 4) + randn(size(n));
Rx = xcorr(xn, 'biased');
Rxx = zeros(1, p+1);
for i = 1 : p+1
Rxx(i) = Rx(N -1 + i); % Rxx存储值为Rxx(0), Rxx(1), ......, Rxx(p) 共p+1个根据样本估计的自相关值
end
a = zeros(p, p);
sigma2 = zeros(1, p);
% 1阶AR模型的初始值
a(1, 1) = - Rxx(2)/Rxx(1);
sigma2(1) = (1 - a(1, 1)^2) * Rxx(1);
% 从p=2阶开始迭代,一直迭代到p阶
for m = 2 : p % m表示当前迭代的AR模型的总的阶数
a(m, m) = -Rxx(m + 1) / sigma2(m - 1); % m阶AR模型总共m个系数中的第m个
for i = 1 : m - 1
a(m, m) = a(m, m) - (a(m-1, i) * Rxx(m +1 - i)) / sigma2(m - 1);
end
for k = 1 : m -1
a(m, k) = a(m - 1, k) + a(m, m) * a(m - 1, m - k)'; % m阶AR模型中的m个系数中的除最后一个系数以外的系数
end
sigma2(m) = (1 - a(m, m) ^ 2) * sigma2(m - 1); % m阶AR模型当前对应的噪声功率
end
[H, w] = freqz(1, [1, a(p, :)], fs);
f = w / (2 * pi) * fs;
H = H * sigma2(m);
plot(f, abs(H).^2);
运行结果:

参考资料:
【1】 张贤达,现代信号处理,第3版
【2】 https://www.bilibili.com/video/BV1wS4y1D7ng/?p=1&vd_source=d7f21d136ca4d445e3fbba3c70c5c40e
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 为DeepSeek添加本地知识库
· 精选4款基于.NET开源、功能强大的通讯调试工具
· DeepSeek智能编程
· 大模型工具KTransformer的安装
· [计算机/硬件/GPU] 显卡