信号处理 | AR模型与Levinson-Durbin递推

模型形式

由高斯白噪声驱动的全极点模型表示如下：
$$\begin{array}{rl}e(n)& =a_{0}x(n)+a_{1}x(n-1)+\cdots +a_{p}x(n-p)\\ & =\sum _{k=0}^p{a}_{k}x(n-k)\end{array}$$
其中：$e(n)$为高斯白噪声，噪声方差为${\sigma }_e^2$，$p$是模型阶数。令$a_0=1$，可以得到如下所示的AR模型：
$$x(n)=-\sum _{k=1}^p{a}_{k}x(n-k)+e(n)$$

$z$域形式为：
$$H(z)=\frac{d_0}{A(z)}=\frac{d_0}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}$$
功率谱密度可表示为：
$$S_{AR}(\omega )=\frac{{\sigma }^2}{|1+{\sum }_{k=1}^p{a}_k{e}^{-j\omega k}{|}^2}$$

模型系数

知道了模型形式，那么如何确定模型系数{$a_k$}？首先给出$x$的自相关函数$r_{xx}$的定义：
$$r_{xx}(i)=E[x(n)x(n-i)]$$
将AR模型带入上式，得到：
$$\begin{array}{rl}{r}_{xx}(i)& =E[(-\sum _{k=1}^p{a}_{k}x(n-k)+e(n))x(n-i)]\\ & =-\sum _{k=1}^p{a}_{k}E[x(n-k)x(n-i)]+E[e(n)x(n-i)]\\ & =\{\begin{array}{ll}-{\sum }_{k=1}^p{a}_k{r}_{xx}(i-k)& i>0\\ -{\sum }_{k=1}^p{a}_k{r}_{xx}(i-k)+{\sigma }_e^2& i=0\end{array}\end{array}$$
即：
$$r_{xx}(i)=\{\begin{array}{ll}-{\sum }_{k=1}^p{a}_k{r}_{xx}(i-k)& i>0\\ -{\sum }_{k=1}^p{a}_k{r}_{xx}(i-k)+{\sigma }_e^2& i=0\end{array}$$
这就是Yule-Walker方程，将上式转换为矩阵相乘形式，其中用到了自相关函数的对称性$r(-i)=r(i)$
$$\left[\begin{array}{cccc}r(0)& r(1)& \cdots & r(p)\\ r(1)& r(0)& \cdots & r(p-1)\\ r(2)& r(1)& \cdots & r(p-2)\\ ⋮& ⋮& & ⋮\\ r(p)& r(p-1)& \cdots & r(1)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ a_2\\ ⋮\\ a_p\end{array}\right]=\left[\begin{array}{c}{\sigma }_e^2\\ 0\\ 0\\ ⋮\\ 0\end{array}\right]$$
一共有$p+1$个方程，未知参数包括$p$个AR模型系数 { a k } \{a_k\} {ak​}以及1个噪声方差${\sigma }_e^2$。观察到上式除了第一个方程，其余方程均与${\sigma }_e^2$无关，因此可以利用其余的$p$个方程求解AR模型系数，将该方程转换为如下形式：
$$\left[\begin{array}{ccc}r(0)& \cdots & r(p-1)\\ r(1)& \cdots & r(p-2)\\ ⋮& & ⋮\\ r(p-1)& \cdots & r(1)\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_p\end{array}\right]=-\left[\begin{array}{c}r(1)\\ r(2)\\ ⋮\\ r(p)\end{array}\right]$$
令$R=\left[\begin{array}{ccc}r(0)& \cdots & r(p-1)\\ r(1)& \cdots & r(p-2)\\ ⋮& & ⋮\\ r(p-1)& \cdots & r(1)\end{array}\right]$，$r=\left[\begin{array}{c}r(1)\\ r(2)\\ ⋮\\ r(p)\end{array}\right]$，$a=\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_p\end{array}\right]$，上式可以表示为$Ra=-r$，所以：
$$a=-R^{-1}r$$

但是当模型阶次较高时，使用矩阵求逆的方法求解模型系数就不合适了，需要更高效的方法。观察到矩阵$R$是一个托普利兹矩阵，可以使用Levinson-Durbin递推来求解AR模型的系数。

Levinson-Durbin递推

已经$x(n)$的$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递推过程如下：

1. 计算$m=1$阶AR模型系数：
   $$\begin{array}{rl}{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{array}$$

2. 递推计算$m=2,3,\cdots ,p$阶AR模型系数：
   $$\begin{array}{rl}{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_m{a}_{m-i}^{(m-1)}\\ {\sigma }_m^2& ={\sigma }_{m-1}^2(1-|k_m{|}^2)\end{array}$$

举一个栗子，假设$r_{xx}(0)=3$，$r_{xx}(1)=2$，$r_{xx}(2)=1$，使用Levinson-Durbin递推2阶AR模型系数：

1. 计算$m=1$阶AR模型系数：
   $$\begin{array}{rl}{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{array}$$

2. 递推计算$m=2$阶AR模型系数：
   $$\begin{array}{rl}{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_2{a}_1^{(1)}=-\frac{4}{5}\\ {\sigma }_2^2& ={\sigma }_1^2(1-|k_2{|}^2)=\frac{8}{5}\end{array}$$

所以2阶AR模型为：
$$x(n)=-\sum _{k=1}^2{a}_{k}x(n-k)=0.8x(n-1)-0.2x(n-2)$$