模型形式
由高斯白噪声驱动的全极点模型表示如下：
e
(
n
)
=
a
0
x
(
n
)
+
a
1
x
(
n
−
1
)
+
⋯
+
a
p
x
(
n
−
p
)
=
∑
k
=
0
p
a
k
x
(
n
−
k
)
\begin{aligned} e(n)&=a_0x(n)+a_1x(n-1)+\cdots+a_px(n-p)\\ &=\sum_{k=0}^pa_kx(n-k) \end{aligned}
e(n)=a0x(n)+a1x(n−1)+⋯+apx(n−p)=k=0∑pakx(n−k)
其中：
e
(
n
)
e(n)
e(n)为高斯白噪声，噪声方差为
σ
e
2
\sigma_e^2
σe2，
p
p
p是模型阶数。令
a
0
=
1
a_0=1
a0=1，可以得到如下所示的AR模型：
x
(
n
)
=
−
∑
k
=
1
p
a
k
x
(
n
−
k
)
+
e
(
n
)
x(n)=-\sum_{k=1}^pa_kx(n-k)+e(n)
x(n)=−k=1∑pakx(n−k)+e(n)
z
z
z域形式为：
H
(
z
)
=
d
0
A
(
z
)
=
d
0
1
+
∑
k
=
1
p
a
k
z
−
k
H(z)=\frac{d_0}{A(z)}=\frac{d_0}{1+\sum_{k=1}^pa_kz^{-k}}
H(z)=A(z)d0=1+∑k=1pakz−kd0
功率谱密度可表示为：
S
A
R
(
ω
)
=
σ
2
∣
1
+
∑
k
=
1
p
a
k
e
−
j
ω
k
∣
2
S_{AR}(\omega)=\frac{\sigma^2}{|1+\sum_{k=1}^pa_ke^{-j\omega k}|^2}
SAR(ω)=∣1+∑k=1pake−jωk∣2σ2
模型系数
知道了模型形式，那么如何确定模型系数{
a
k
a_k
ak}？首先给出
x
x
x的自相关函数
r
x
x
r_{xx}
rxx的定义：
r
x
x
(
i
)
=
E
[
x
(
n
)
x
(
n
−
i
)
]
r_{xx}(i)=E[x(n)x(n-i)]
rxx(i)=E[x(n)x(n−i)]
将AR模型带入上式，得到：
r
x
x
(
i
)
=
E
[
(
−
∑
k
=
1
p
a
k
x
(
n
−
k
)
+
e
(
n
)
)
x
(
n
−
i
)
]
=
−
∑
k
=
1
p
a
k
E
[
x
(
n
−
k
)
x
(
n
−
i
)
]
+
E
[
e
(
n
)
x
(
n
−
i
)
]
=
{
−
∑
k
=
1
p
a
k
r
x
x
(
i
−
k
)
i
>
0
−
∑
k
=
1
p
a
k
r
x
x
(
i
−
k
)
+
σ
e
2
i
=
0
\begin{aligned} r_{xx}(i)&=E[(-\sum_{k=1}^pa_kx(n-k)+e(n))x(n-i)]\\ &=-\sum_{k=1}^pa_kE[x(n-k)x(n-i)]+E[e(n)x(n-i)]\\ &=\begin{cases}-\sum_{k=1}^pa_kr_{xx}(i-k)&i>0\\ -\sum_{k=1}^pa_kr_{xx}(i-k)+\sigma_e^2&i=0 \end{cases} \end{aligned}
rxx(i)=E[(−k=1∑pakx(n−k)+e(n))x(n−i)]=−k=1∑pakE[x(n−k)x(n−i)]+E[e(n)x(n−i)]={−∑k=1pakrxx(i−k)−∑k=1pakrxx(i−k)+σe2i>0i=0
即：
r
x
x
(
i
)
=
{
−
∑
k
=
1
p
a
k
r
x
x
(
i
−
k
)
i
>
0
−
∑
k
=
1
p
a
k
r
x
x
(
i
−
k
)
+
σ
e
2
i
=
0
r_{xx}(i)=\begin{cases}-\sum_{k=1}^pa_kr_{xx}(i-k)&i>0\\ -\sum_{k=1}^pa_kr_{xx}(i-k)+\sigma_e^2&i=0\end{cases}
rxx(i)={−∑k=1pakrxx(i−k)−∑k=1pakrxx(i−k)+σe2i>0i=0
这就是Yule-Walker方程，将上式转换为矩阵相乘形式，其中用到了自相关函数的对称性
r
(
−
i
)
=
r
(
i
)
r(-i)=r(i)
r(−i)=r(i)
[
r
(
0
)
r
(
1
)
⋯
r
(
p
)
r
(
1
)
r
(
0
)
⋯
r
(
p
−
1
)
r
(
2
)
r
(
1
)
⋯
r
(
p
−
2
)
⋮
⋮
⋮
r
(
p
)
r
(
p
−
1
)
⋯
r
(
1
)
]
[
1
a
1
a
2
⋮
a
p
]
=
[
σ
e
2
0
0
⋮
0
]
\begin{bmatrix} r(0)&r(1)&\cdots&r(p)\\ r(1)&r(0)&\cdots&r(p-1)\\ r(2)&r(1)&\cdots&r(p-2)\\ \vdots&\vdots&&\vdots\\ r(p)&r(p-1)&\cdots&r(1) \end{bmatrix} \begin{bmatrix} 1\\a_1\\a_2\\\vdots\\a_p \end{bmatrix} =\begin{bmatrix}\sigma^2_e\\0\\0\\\vdots\\0\end{bmatrix}
⎣⎢⎢⎢⎢⎢⎡r(0)r(1)r(2)⋮r(p)r(1)r(0)r(1)⋮r(p−1)⋯⋯⋯⋯r(p)r(p−1)r(p−2)⋮r(1)⎦⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡1a1a2⋮ap⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡σe200⋮0⎦⎥⎥⎥⎥⎥⎤
一共有
p
+
1
p+1
p+1个方程，未知参数包括
p
p
p个AR模型系数
{
a
k
}
\{a_k\}
{ak}以及1个噪声方差
σ
e
2
\sigma_e^2
σe2。观察到上式除了第一个方程，其余方程均与
σ
e
2
\sigma_e^2
σe2无关，因此可以利用其余的
p
p
p个方程求解AR模型系数，将该方程转换为如下形式：
[
r
(
0
)
⋯
r
(
p
−
1
)
r
(
1
)
⋯
r
(
p
−
2
)
⋮
⋮
r
(
p
−
1
)
⋯
r
(
1
)
]
[
a
1
a
2
⋮
a
p
]
=
−
[
r
(
1
)
r
(
2
)
⋮
r
(
p
)
]
\begin{bmatrix} r(0)&\cdots&r(p-1)\\ r(1)&\cdots&r(p-2)\\ \vdots&&\vdots\\ r(p-1)&\cdots&r(1) \end{bmatrix} \begin{bmatrix} a_1\\a_2\\\vdots\\a_p \end{bmatrix} =-\begin{bmatrix}r(1)\\r(2)\\\vdots\\r(p)\end{bmatrix}
⎣⎢⎢⎢⎡r(0)r(1)⋮r(p−1)⋯⋯⋯r(p−1)r(p−2)⋮r(1)⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡a1a2⋮ap⎦⎥⎥⎥⎤=−⎣⎢⎢⎢⎡r(1)r(2)⋮r(p)⎦⎥⎥⎥⎤
令
R
=
[
r
(
0
)
⋯
r
(
p
−
1
)
r
(
1
)
⋯
r
(
p
−
2
)
⋮
⋮
r
(
p
−
1
)
⋯
r
(
1
)
]
R=\begin{bmatrix} r(0)&\cdots&r(p-1)\\ r(1)&\cdots&r(p-2)\\ \vdots&&\vdots\\ r(p-1)&\cdots&r(1) \end{bmatrix}
R=⎣⎢⎢⎢⎡r(0)r(1)⋮r(p−1)⋯⋯⋯r(p−1)r(p−2)⋮r(1)⎦⎥⎥⎥⎤，
r
=
[
r
(
1
)
r
(
2
)
⋮
r
(
p
)
]
r=\begin{bmatrix}r(1)\\r(2)\\\vdots\\r(p)\end{bmatrix}
r=⎣⎢⎢⎢⎡r(1)r(2)⋮r(p)⎦⎥⎥⎥⎤，
a
=
[
a
1
a
2
⋮
a
p
]
a=\begin{bmatrix} a_1\\a_2\\\vdots\\a_p \end{bmatrix}
a=⎣⎢⎢⎢⎡a1a2⋮ap⎦⎥⎥⎥⎤，上式可以表示为
R
a
=
−
r
Ra=-r
Ra=−r，所以：
a
=
−
R
−
1
r
a=-R^{-1}r
a=−R−1r
但是当模型阶次较高时，使用矩阵求逆的方法求解模型系数就不合适了，需要更高效的方法。观察到矩阵
R
R
R是一个托普利兹矩阵，可以使用Levinson-Durbin递推来求解AR模型的系数。
Levinson-Durbin递推
已经
x
(
n
)
x(n)
x(n)的
p
+
1
p+1
p+1个自相关函数
{
r
x
x
(
0
)
,
r
x
x
(
1
)
,
⋯
,
r
x
x
(
p
)
}
\{r_{xx}(0),r_{xx}(1),\cdots,r_{xx}(p)\}
{rxx(0),rxx(1),⋯,rxx(p)}，Levinson-Durbin递推过程如下：
-
计算
m
=
1
m=1
m=1阶AR模型系数：
a
1
(
1
)
=
−
r
x
x
(
1
)
r
x
x
(
0
)
σ
1
2
=
r
x
x
(
0
)
−
∣
r
x
x
(
1
)
∣
2
r
x
x
(
0
)
\begin{aligned} a_1^{(1)}&=-\frac{r_{xx}(1)}{r_{xx}(0)}\\ \sigma_1^2&=r_{xx}(0)-\frac{|r_{xx}(1)|^2}{r_{xx}(0)} \end{aligned}
a1(1)σ12=−rxx(0)rxx(1)=rxx(0)−rxx(0)∣rxx(1)∣2
-
递推计算
m
=
2
,
3
,
⋯
,
p
m=2,3,\cdots,p
m=2,3,⋯,p阶AR模型系数：
k
m
=
r
x
x
(
m
)
+
∑
i
=
1
m
−
1
a
i
(
m
−
1
)
r
x
x
(
m
−
i
)
σ
m
−
1
2
a
m
(
m
)
=
k
m
a
i
(
m
)
=
a
i
(
m
−
1
)
+
k
m
a
m
−
i
(
m
−
1
)
σ
m
2
=
σ
m
−
1
2
(
1
−
∣
k
m
∣
2
)
\begin{aligned} k_m&=\frac{r_{xx}(m)+\sum_{i=1}^{m-1}a_i^{(m-1)}r_{xx}(m-i)}{\sigma_{m-1}^2}\\ a_m^{(m)}&=k_m\\ a_i^{(m)}&=a_i^{(m-1)}+k_ma_{m-i}^{(m-1)}\\ \sigma_m^2&=\sigma_{m-1}^2(1-|k_m|^2)\\ \end{aligned}
kmam(m)ai(m)σm2=σm−12rxx(m)+∑i=1m−1ai(m−1)rxx(m−i)=km=ai(m−1)+kmam−i(m−1)=σm−12(1−∣km∣2)
举一个栗子，假设
r
x
x
(
0
)
=
3
r_{xx}(0)=3
rxx(0)=3，
r
x
x
(
1
)
=
2
r_{xx}(1)=2
rxx(1)=2，
r
x
x
(
2
)
=
1
r_{xx}(2)=1
rxx(2)=1，使用Levinson-Durbin递推2阶AR模型系数：
-
计算
m
=
1
m=1
m=1阶AR模型系数：
a
1
(
1
)
=
−
r
x
x
(
1
)
r
x
x
(
0
)
=
−
2
3
σ
1
2
=
r
x
x
(
0
)
−
∣
r
x
x
(
1
)
∣
2
r
x
x
(
0
)
=
5
3
\begin{aligned} a_1^{(1)}&=-\frac{r_{xx}(1)}{r_{xx}(0)}=-\frac{2}{3}\\ \sigma_1^2&=r_{xx}(0)-\frac{|r_{xx}(1)|^2}{r_{xx}(0)}=\frac{5}{3} \end{aligned}
a1(1)σ12=−rxx(0)rxx(1)=−32=rxx(0)−rxx(0)∣rxx(1)∣2=35
-
递推计算
m
=
2
m=2
m=2阶AR模型系数：
k
2
=
r
x
x
(
2
)
+
a
1
(
1
)
r
x
x
(
1
)
σ
1
2
=
1
5
a
2
(
2
)
=
k
2
=
1
5
a
1
(
2
)
=
a
1
(
1
)
+
k
2
a
1
(
1
)
=
−
4
5
σ
2
2
=
σ
1
2
(
1
−
∣
k
2
∣
2
)
=
8
5
\begin{aligned} k_2&=\frac{r_{xx}(2)+a_1^{(1)}r_{xx}(1)}{\sigma_{1}^2}=\frac{1}{5}\\ a_2^{(2)}&=k_2=\frac{1}{5}\\ a_1^{(2)}&=a_1^{(1)}+k_2a_{1}^{(1)}=-\frac{4}{5}\\ \sigma_2^2&=\sigma_{1}^2(1-|k_2|^2)=\frac{8}{5}\\ \end{aligned}
k2a2(2)a1(2)σ22=σ12rxx(2)+a1(1)rxx(1)=51=k2=51=a1(1)+k2a1(1)=−54=σ12(1−∣k2∣2)=58
所以2阶AR模型为：
x
(
n
)
=
−
∑
k
=
1
2
a
k
x
(
n
−
k
)
=
0.8
x
(
n
−
1
)
−
0.2
x
(
n
−
2
)
x(n)=-\sum_{k=1}^2a_kx(n-k)=0.8x(n-1)-0.2x(n-2)
x(n)=−k=1∑2akx(n−k)=0.8x(n−1)−0.2x(n−2)
本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)
