The Levinson Recursion_statistical digital signal processing and modeling-CSDN博客

本文链接：https://blog.csdn.net/qq_39599295/article/details/108935601

Reference:
Slides of EE4C03, TUD
Hayes M H. Statistical digital signal processing and modeling

Content

Development of the Recursion

All-pole modeling using Prony’s method or the autocorrelation method requires that we solve the normal equations which, for a $p$ th-order model, are
$r_x(k)+\sum_{l=1}^pa_p(l)r_x(k-l)=0;\quad k=1,2,\cdots ,p\tag{D.1}$
where the modeling error is
$\epsilon_p=r_x(0)+\sum_{l=l}^pa_p(l)r_x(l)\tag{D.2}$
Combining $(D . 1)$ and $(D . 2)$ into matrix form we have
$\left[\begin{array}{ccccc} r_{x}(0) & r_{x}^{*}(1) & r_{x}^{*}(2) & \cdots & r_{x}^{*}(p) \\ r_{x}(1) & r_{x}(0) & r_{x}^{*}(1) & \cdots & r_{x}^{*}(p-1) \\ r_{x}(2) & r_{x}(1) & r_{x}(0) & \cdots & r_{x}^{*}(p-2) \\ \vdots & \vdots & \vdots & & \vdots \\ r_{x}(p) & r_{x}(p-1) & r_{x}(p-2) & \cdots & r_{x}(0) \end{array}\right]\left[\begin{array}{c} 1 \\ a_{p}(1) \\ a_{p}(2) \\ \vdots \\ a_{p}(p) \end{array}\right]=\epsilon_{p}\left[\begin{array}{c} 1 \\ 0 \\ 0 \\ \vdots \\ 0 \end{array}\right] \tag{D.3}$
which is a set of $p + 1$ linear equations in the $p + 1$ unknowns $a_p(1),a_p(2),\cdots, a_p(p)$ and $\epsilon_p$ . Equivalently, $(D . 3)$ may be written as
$\mathbf R_p \mathbf a_p=\epsilon_p \mathbf u_1\tag{D.4}$

The Levinson-Durbin recursion for solving (D.4) is an algorithm that is recursive in the model order. In other words, the coefficients of the $(j + 1)$ st-order all-pole model, $\mathbf a_{j+1}$ , are found from the coefficients of the $j$ -pole model, $\mathbf a_j$ .

Let $a_j (i)$ be the solution to the $j$ th-order normal equations
$\left[\begin{array}{ccccc} r_{x}(0) & r_{x}^{*}(1) & r_{x}^{*}(2) & \cdots & r_{x}^{*}(j) \\ r_{x}(1) & r_{x}(0) & r_{x}^{*}(1) & \cdots & r_{x}^{*}(j-1) \\ r_{x}(2) & r_{x}(1) & r_{x}(0) & \cdots & r_{x}^{*}(j-2) \\ \vdots & \vdots & \vdots & & \vdots \\ r_{x}(j) & r_{x}(j-1) & r_{x}(j-2) & \cdots & r_{x}(0) \end{array}\right]\left[\begin{array}{c} 1 \\ a_{p}(1) \\ a_{p}(2) \\ \vdots \\ a_{p}(j) \end{array}\right]=\left[\begin{array}{c} \epsilon_{j} \\ 0 \\ 0 \\ \vdots \\ 0 \end{array}\right] \tag{D.5}$
which, in matrix notation is
$\mathbf R_j \mathbf a_j=\epsilon_j \mathbf u_1\tag{D.6}$
Given $\mathbf a_j$ , we want to derive the solution to the $(j + 1)$ st-order normal equations,
$\mathbf R_{j+1} \mathbf a_{j+1}=\epsilon_{j+1} \mathbf u_1\tag{D.7}$
Suppose we append a zero to the vector $\mathbf a_j$ and multiply the resulting vector by $\mathbf R_{j+1}$ . The result is
$\left[\begin{array}{ccccc} r_{x}(0) & r_{x}^{*}(1) & r_{x}^{*}(2) & \cdots & r_{x}^{*}(j) & r_{x}^{*}(j+1) \\ r_{x}(1) & r_{x}(0) & r_{x}^{*}(1) & \cdots & r_{x}^{*}(j-1) & r_{x}^{*}(j) \\ r_{x}(2) & r_{x}(1) & r_{x}(0) & \cdots & r_{x}^{*}(j-2) & r_{x}^{*}(j-1) \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ r_{x}(j) & r_{x}(j-1) & r_{x}(j-2) & \cdots & r_{x}(0) & r_{x}^{*}(1) \\ r_{x}(j+1) & r_{x}(j) & r_{x}(j-1) & \cdots & r_{x}(1) & r_{x}(0) \end{array}\right]\left[\begin{array}{c} 1 \\ a_{j}(1) \\ a_{j}(2) \\ \vdots \\ a_{j}(j) \\ 0 \end{array}\right]=\left[\begin{array}{c} \epsilon_{j} \\ 0 \\ 0 \\ \vdots \\ 0 \\ \gamma_{j} \end{array}\right] \tag{D.8}$
where the parameter $\gamma_j$ is
$\gamma_j=r_x(j+1)+\sum_{i=1}^j a_j(i)r_x(j+1-i)\tag{D.9}$
Note that if $\gamma_{j}=0,$ then the right side of $(D . 8)$ is a scaled unit vector and $\mathbf{a}_{j+1}=$ $\left[1, a_{j}(1), \ldots, a_{j}(j), 0\right]^{T}$ is the solution to the $(j + 1)$ st-order normal equations $(D . 7) .$ In general, however, $\gamma_{j} \neq 0$ and $\left[1, a_{j}(1), \ldots, a_{j}(j), 0\right]^{T}$ is not the solution to $(D . 7)$ .

The key step in the derivation of the Levinson-Durbin recursion is to note that the Hermitian Toeplitz property of $R_{j+1}$ allows us to rewrite $(D . 8)$ in the equivalent form
$\left[\begin{array}{ccccc} r_{x}(0) & r_{x}(1) & r_{x}(2) & \cdots & r_{x}(j) & r_{x}(j+1) \\ r_{x}^{*}(1) & r_{x}(0) & r_{x}(1) & \cdots & r_{x}(j-1) & r_{x}(j) \\ r_{x}^{*}(2) & r_{x}^{*}(1) & r_{x}(0) & \cdots & r_{x}(j-2) & r_{x}(j-1) \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ r_{x}^{*}(j) & r_{x}^{*}(j-1) & r_{x}^{*}(j-2) & \cdots & r_{x}(0) & r_{x}(1) \\ r_{x}^{*}(j+1) & r_{x}^{*}(j) & r_{x}^{*}(j-1) & \cdots & r_{x}^{*}(1) & r_{x}(0) \end{array}\right]\left[\begin{array}{c} 0 \\ a_{j}(j) \\ a_{j}(j-1) \\ \vdots \\ a_{j}(1) \\ 1 \end{array}\right]=\left[\begin{array}{c} \gamma_{j} \\ 0 \\ 0 \\ \vdots \\ 0 \\ \epsilon_{j} \end{array}\right] \tag{D.10}$
Taking the complex conjugate of $(D . 10)$
$\left[\begin{array}{ccccc} r_{x}(0) & r_{x}^{*}(1) & r_{x}^{*}(2) & \cdots & r_{x}^{*}(j) & r_{x}^{*}(j+1) \\ r_{x}(1) & r_{x}(0) & r_{x}^{*}(1) & \cdots & r_{x}^{*}(j-1) & r_{x}^{*}(j) \\ r_{x}(2) & r_{x}(1) & r_{x}(0) & \cdots & r_{x}^{*}(j-2) & r_{x}^{*}(j-1) \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ r_{x}(j) & r_{x}(j-1) & r_{x}(j-2) & \cdots & r_{x}(0) & r_{x}^{*}(1) \\ r_{x}(j+1) & r_{x}(j) & r_{x}(j-1) & \cdots & r_{x}(1) & r_{x}(0) \end{array}\right]\left[\begin{array}{c} 0 \\ a_{j}^*(j) \\ a_{j}^*(j-1) \\ \vdots \\ a_{j}^*(1) \\ 1 \end{array}\right]=\left[\begin{array}{c} \gamma_{j}^* \\ 0 \\ 0 \\ \vdots \\ 0 \\ \epsilon_{j}^* \end{array}\right] \tag{D.11}$
Although $\left[1, a_{j}(1), \ldots, a_{j}(j), 0\right]^{T}$ is not the solution to $(D . 7)$ , the linear combination of $\left[1, a_{j}(1), \ldots, a_{j}(j), 0\right]^{T}$ and $\left[0, a_{j}^*(j), \ldots, a_{j}^*(1), 1\right]^{T}$ can be the solution to $(D . 7)$ . To see this,
$\mathbf{R}_{j+1}\left\{\left[\begin{array}{c} 1 \\ a_{j}(1) \\ a_{j}(2) \\ \vdots \\ a_{j}(j) \\ 0 \end{array}\right]+\Gamma_{j+1}\left[\begin{array}{c} 0 \\ a_{j}^{*}(j) \\ a_{j}^{*}(j-1) \\ \vdots \\ a_{j}^{*}(1) \\ 1 \end{array}\right]\right\}=\left[\begin{array}{c} \epsilon_{j} \\ 0 \\ 0 \\ \vdots \\ 0 \\ \gamma_{j} \end{array}\right]+\Gamma_{j+1}\left[\begin{array}{c} \gamma_{j}^{*} \\ 0 \\ 0 \\ \vdots \\ 0 \\ \epsilon_{j}^{*} \end{array}\right]\tag{D.12}$
If we set
$\Gamma_{j+1}=-\frac{\gamma_j}{\epsilon_j^*} \tag{D.13}$
then $(D . 12)$ becomes
$\mathbf R_{j+1}\mathbf a_{j+1}=\epsilon_{j+1}\mathbf u_1$
where
$\mathbf{a}_{j+1}=\left[\begin{array}{c} 1 \\ a_{j}(1) \\ a_{j}(2) \\ \vdots \\ a_{j}(j) \\ 0 \end{array}\right]+\Gamma_{j+1}\left[\begin{array}{c} 0 \\ a_{j}^{*}(j) \\ a_{j}^{*}(j-1) \\ \vdots \\ a_{j}^{*}(1) \\ 1 \end{array}\right]\tag{D.14}$
Furthermore,
$\epsilon_{j+1}=\epsilon_j+\Gamma_{j+1}\gamma_j^*=\epsilon_j[1-|\Gamma_{j+1}|^2]\tag{D.15}$
The whole process can also be presented as
$\mathbf{R}_{j+1} \left[\begin{array}{ccc} 1 & 0 \\ a_{j}(1) & a_{j}^*(j) \\ \vdots & \vdots \\ a_{j}(j) & a_{j}^*(1) \\ 0 & 1 \end{array}\right] \left[\begin{array}{ccc} 1 & \Gamma_{j+1}^* \\ \Gamma_{j+1} & 1 \end{array}\right]=\mathbf{R}_{j+1} \left[\begin{array}{ccc} 1 & a_{j+1}^*(j+1) \\ a_{j+1}(1) & a_{j+1}^*(j) \\ \vdots & \vdots \\ a_{j+1}(j) & a_{j+1}^*(1) \\ a_{j+1}(j+1) & 1 \end{array}\right] =\left[\begin{array}{cc} \epsilon_{j} & \gamma_{j}^* \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \\ \gamma_{j} & \epsilon_{j}^* \end{array}\right]\left[\begin{array}{ccc} 1 & \Gamma_{j+1}^* \\ \Gamma_{j+1} & 1 \end{array}\right]=\left[\begin{array}{cc} \epsilon_{j+1} & 0\\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \\ 0 & \epsilon_{j+1}^* \end{array}\right]\tag{D.16}$
All that is required to complete the recursion is to define the conditions necessary to initialize the recursion:
$a_0(0)=1,\quad \epsilon_0=r_x(0) \tag{D.17}$
In summary, the steps of the Levinson-Durbin recursion are as follows:
$\left\{r_{x}(0), r_{x}(1), \ldots, r_{x}(p)\right\} \stackrel{L E V}{\longrightarrow}\left\{\begin{array}{c} \Gamma_{1}, \Gamma_{2}, \ldots, \Gamma_{p}, \epsilon_{p}\\ a_{p}(1), a_{p}(2), \ldots, a_{p}(p), b(0) \end{array}\right.$

在这里插入图片描述

Alternatively, using the notation in $(D . 16)$ ,

Get $\epsilon_j,\gamma_j$ from
$\mathbf{R}_{j+1} \left[\begin{array}{ccc} 1 & 0 \\ a_{j}(1) & a_{j}^*(j) \\ \vdots & \vdots \\ a_{j}(j) & a_{j}^*(1) \\ 0 & 1 \end{array}\right] =\left[\begin{array}{cc} \epsilon_{j} & \gamma_{j}^* \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \\ \gamma_{j} & \epsilon_{j}^* \end{array}\right]$
Obtain $\Gamma_{j+1}=-\gamma_j/\epsilon_j$ .
Obtain $\mathbf a_{j+1}$ from
$\left[\begin{array}{ccc} 1 & 0 \\ a_{j}(1) & a_{j}^*(j) \\ \vdots & \vdots \\ a_{j}(j) & a_{j}^*(1) \\ 0 & 1 \end{array}\right] \left[\begin{array}{ccc} 1 & \Gamma_{j+1}^* \\ \Gamma_{j+1} & 1 \end{array}\right]= \left[\begin{array}{ccc} 1 & a_{j+1}^*(j+1) \\ a_{j+1}(1) & a_{j+1}^*(j) \\ \vdots & \vdots \\ a_{j+1}(j) & a_{j+1}^*(1) \\ a_{j+1}(j+1) & 1 \end{array}\right]$

The Lattice Filter

Define the reciprocal vector $\mathbf a_j^R$ , which is the vector that is formed by reversing the order of the elements in $\mathbf a_j$ and taking the complex conjugate,
$\mathbf{a}_{j}=\left[\begin{array}{c} 1 \\ a_{j}(1) \\ a_{j}(2) \\ \vdots \\ a_{j}(j-1) \\ a_{j}(j) \end{array}\right] \Longrightarrow\left[\begin{array}{c} a_{j}^{*}(j) \\ a_{j}^{*}(j-1) \\ a_{j}^{*}(j-2) \\ \vdots \\ a_{j}^{*}(1) \\ 1 \end{array}\right]=\mathbf{a}_{j}^{R} \tag{LF.1}$
or
$a_j^R(i)=a_j^*(j-i),\text { for }i=0,1,\cdots,j \tag{LF.2}$

FIR lattice filter

From $(D . 16)$ , we have
$\left[\begin{array}{ccc} 1 & 0 \\ a_{j}(1) & a_{j}^*(j) \\ \vdots & \vdots \\ a_j(i) & a_j^*(j-i+1) \\ \vdots & \vdots \\ a_{j}(j) & a_{j}^*(1) \\ 0 & 1 \end{array}\right] \left[\begin{array}{ccc} 1 & \Gamma_{j+1}^* \\ \Gamma_{j+1} & 1 \end{array}\right]= \left[\begin{array}{ccc} 1 & a_{j+1}^*(j+1) \\ a_{j+1}(1) & a_{j+1}^*(j) \\ \vdots & \vdots \\ a_{j+1}(i) & a_{j+1}^*(j+1-i) \\ \vdots & \vdots \\ a_{j+1}(j) & a_{j+1}^*(1) \\ a_{j+1}(j+1) & 1 \end{array}\right]$

i.e.,
$\begin{aligned} & a_{j+1}(n)=a_j(n)+\Gamma_{j+1}a_j^R(n-1)\\ & a_{j+1}^R(n)=a_j^R(n-1)+\Gamma^*_{j+1}a_j(n) \end{aligned} \tag{LF.3}$

在这里插入图片描述

We can define filter functions corresponding to the filter coefficients
$\left\{\begin{array}{l} A_{j}(z)=1+a_{j}(1) z^{-1}+\cdots+a_{j}(j) z^{-j} \\ A_{j}^{R}(z)=a_{j}^*(j)+a_{j}^*(j-1) z^{-1}+\cdots+z^{-j} \end{array}\right. \tag{LF.4}$
which are the $z$ -transform of sequences $\mathbf a_j=[1,a_j(1),\cdots,a_j(j)]^T$ and $\mathbf a_j^R=[a_j^*(j),a_j^*(j-1),\cdots,1]^T$ .

According to $(L F . 2)$ , it can be verified that
$A_j^R(z)=z^{-j}A_j^*(1/z^*) \tag{LF.5}$
From $(L F . 3)$ , we obtain
$\begin{aligned} & A_{j+1}(z)=A_j(z)+z^{-1}\Gamma_{j+1}A_j^R(z)\\ & A_{j+1}^R(z)=z^{-1}A_j^R(z)+\Gamma_{j+1}^*A_j(z) \end{aligned} \tag{LF.6}$
(Note that the notations of the slides is different from the book. $\rho$ equals to $-\Gamma$ and the slides only consider real valued signals)

在这里插入图片描述

Thus, the resulting system has impulse response $A_p(z)$ (and also $A_P^R(z)$ , not used).

If we use this system with input $x [n]$ , we obtain the prediction error sequence $v_P[n]$ : (Note that in this process, we only need to know $\rho_i$ , $a_p(n)$ is not necessary)

在这里插入图片描述

The filter structure is known as a FIR lattice filter: the filter response is $A_p(z)$ .

IIR lattice filter

Rearrange $(L F . 6)$ as
$\begin{aligned} & A_{j}(z)=A_{j+1}(z)-z^{-1}\Gamma_{j+1}A_j^R(z)\\ & A_{j+1}^R(z)=z^{-1}A_j^R(z)+\Gamma_{j+1}^*A_j(z) \end{aligned} \tag{LF.7}$

在这里插入图片描述

This allows to compute $x [n]$ from $v_p[n]$ : the filter response is $\frac{1}{A_p(z)}$ . This is an IIR filter and the filter coefficients of $\frac{1}{A_p(z)}$ are not explicitly computed.

From the correlation sequence $\{r_x(0),\cdots ,r_x(p)\}$ of $x [n]$ , we can obtain $\rho_1,\cdots,\rho_p$ and the variance $\epsilon_p$ of $v_p[n]$ . If we replace $v_p[n]$ by any white noise sequence with variance $\epsilon_p$ , then the resulting output signal is a random process with the same correlation sequence $\{r_x(0),\cdots ,r_x(p)\}$ as $x [n]$ .

Properties

more in book Chapter 5.2.3

To show that the Levinson recursion works and does not break down, we need to prove that we always can compute a suitable $\rho$ . This will follow from the main assumption:
$\mathbf{R}_{x}\succ0 \text{ for all }p.$

$\epsilon_{p}>0 .$

This follows from the YW equations:
$\left[\begin{array}{c} 1 \\ a_{p}(1) \\ a_{p}(2) \\ \vdots \\ a_{p}(p) \end{array}\right]=\left[\begin{array}{ccccc} r_{x}(0) & r_{x}^*(1) & r_{x}^*(2) & \cdots & r_{x}^*(p) \\ r_{x}(1) & r_{x}(0) & r_{x}^*(1) & \cdots & r_{x}^*(p-1) \\ r_{x}(2) & r_{x}(1) & r_{x}(0) & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & r_{x}^*(1) \\ r_{x}(p) & r_{x}(p-1) & \cdots & r_{x}(1) & r_{x}(0) \end{array}\right]^{-1}\left[\begin{array}{c} \epsilon_{p} \\ 0 \\ 0 \\ \vdots \\ 0 \end{array}\right]$
In particular, $1=\left[\mathrm{R}_{x}^{-1}\right]_{0,0} \epsilon_{p} .$ since $\mathrm{R}_{x}$ is strictly positive definite (hence also invertible), the inverse $\mathrm{R}_{x}^{-1}$ exists and is also strictly positive definite. This implies that $\left[\mathbf{R}_{x}^{-1}\right]_{0,0}>0 .$ Hence $\epsilon_{p}>0$ for any $p$ . In particular, $\epsilon_{p} \neq 0,$ so that we can compute $\rho_{\rho+1}$ .
$\left|\rho_{p+1}\right|<1$
The update equation is
$\epsilon_{p+1}=\epsilon_{p}-\rho_{p+1} \gamma_{p}=\epsilon_{p}-\frac{\gamma_{p}^{2}}{\epsilon_{p}}=\frac{\epsilon_{p}^{2}-\gamma_{p}^{2}}{\epsilon_{p}}>0 \Rightarrow\left|\epsilon_{p}\right|>\left|\gamma_{p}\right|$
so that $\left|\rho_{p+1}\right|<1$ . From this we also see that $\epsilon_{p+1} \leq \epsilon_{p}:$ the modeling/prediction error always decreases.

Special case:

If we have a true AR process of order $p,$ then we will find $\epsilon_{p+1}=\epsilon_{p}$ and $\gamma_{p}=0 .$ At this point the recursion stops (we can take $\rho_{p+1}=0,$ so that the AR coefficients do not change anymore).

Special case:

If $\left|\rho_{p+1}\right|=1,$ then $\epsilon_{p+1}=0 .$ The prediction error is zero, $x [n]$ can be exactly predicted from its past, the process is called deterministic. This can occur only if $\mathbf{R}_{x}$ is singular (thus the matrix is not strictly positive definite).
The roots of $A_P(z)$ lie inside the unit circle, i.e., the IIR filter is stable.

The Schur Algorithm

The Levinson algorithm involves two times $2 p$ multiplications and $p$ additions. With $p$ computational processors, this work can be done in parallel. However, the computation of $\gamma_{p}(p$ additions) is not easily parallelized.
The Schur algorithm is an alternative to Levinson to solve the same equations. It is based on the idea that we do not need the filter coefficients $\left\{a_{p}(1), \cdots a_{p}(p)\right\}$ for all $p,$ the reflection coefficients $\left\{\rho_{1}, \cdots, \rho_{p}\right\}$ completely specify the filter.

The development of the Schur recursion begins with the autocorrelation normal equations for a $j$ th-order model (see $(D . 1)$ )
$r_x(k)+\sum_{l=1}^ja_j(l)r_x(k-l)=0;\quad k=1,2,\cdots ,j$
If we set $a_j(0)=1$ then the $j$ th-order normal equations may be written as
$\sum_{l=0}^ja_j(l)r_x(k-l)=0;\quad k=1,2,\cdots ,j$
which we will refer to as the orthogonality condition.

By introducing new variables $g_j(k)$ and $g_j^R(k)$ , we can avoid the computation of the filter coefficients $a_j(k)$ .

Define
$g_j(k)=\sum_{l=0}^j a_j(l)r_x(k-l)=a_j(k)*r_x(k) \tag{S.1}$
i.e., $g_j(k)$ is the response of the filter $A_j(z)$ to the input $r_x(k)$ . Then the orthogonality condition can be presented as
$g_j(k)=0;\quad k=1,2,\cdots,j\tag{S.2}$
In addition, from $(D . 2)$ , we obtain that $g_j(0)$ is equal to the $j$ th-order modeling error:
$g_j(0)=\sum_{l=0}^ja_j(l)r_x(l)=\epsilon_j \tag{S.3}$
(Note that the length of $\mathbf g_j$ is not $j + 1$ but $p + 1$ .)

Similarly, we can define
$g_{j}^{R}(k)=\sum_{l=0}^{j} a_{j}^{R}(l) r_{x}(k-l)=a_{j}^{R}(k) * r_{x}(k)\tag{S.4}$
Thus, $g_{j}^{R}(k)$ is the response of the filter $A_{j}^{R}(z)$ to the input $r_{x}(k) .$ since $a_{j}^{R}(k)=a_{j}^{*}(j-k)$ then
$g_{j}^{R}(k)=\sum_{l=0}^{j} a_{j}^{*}(j-l) r_{x}(k-l)=\sum_{l=0}^{j} a_{j}^{*}(l) r_{x}(k-j+l)\tag{S.5}$
Using the conjugate symmetry of the autocorrelation sequence, $(S . 5)$ becomes
$g_{j}^{R}(k)=\sum_{l=0}^{j} a_{j}^{*}(l) r_{x}^{*}([j-k]-l)=g_{j}^{*}(j-k)\tag{S.6}$
Therefore, it follows from $(S . 2)$ and $(S . 3)$ that
$g_j^R(k)=0;\quad k=0,1,\cdots,j-1 \tag{S.7}$
and
$g_j^R(j)=\epsilon_j \tag{S.8}$
To see this in matrix form,
$\left[\begin{matrix}\mathbf g_j^T\\(\mathbf g_j^R)^T\end{matrix}\right]=\left[\begin{matrix}\epsilon_j & \cdots & 0 & 0& g_j(j+1)& \cdots& g_j(p)\\0 & \cdots & 0 & \epsilon_j& g_j^R(j+1)& \cdots& g^R_j(p)\end{matrix}\right]\tag{S.9}$

The next step is to use the Levinson-Durbin recursion to show how the sequences $g_j(k)$ and $g_j^R(k)$ may be updated to form $g_{j+1}(k)$ and $g_{j+1}^R (k)$ .

Using the Levinson order-update equation $(L F . 3)$ , we have
$\begin{aligned} &g_{j+1}(k)=a_{j+1}(k)*r_x(k)=[a_j(k)+\Gamma_{j+1}a_j^R(k-1)]*r_x(k)=g_j(k)+\Gamma_{j+1}g_j^R(k-1)\\ &g^R_{j+1}(k)=a^R_{j+1}(k)*r_x(k)=[a_j^R(k-1)+\Gamma_{j+1}^*a_j(k)]*r_x(k)=g_j^R(k-1)+\Gamma_{j+1}^*g_j(k) \end{aligned} \tag{S.10}$
Note that the recursive equations for $g_j(k),g^R_j(k)$ $(S . 10)$ and the recursive equations for $a_j(k),a^R_j(k)$ $(L F . 3)$ are identical. The only difference between the two recursions is in the initial condition. Specifically, for $a_j(k)$ we have $a_0(k)=a_0^R(k)=\delta(k)$ , whereas for $g_j(k)$ we have $g_0(k)=g_0^R(k)=r_x(k)$ .

Taking $z$ -transform of $(S . 10)$ and putting them in matrix form gives
$\left[\begin{matrix}G_{j+1}(z)\\ G_{j+1}^R(z)\end{matrix}\right]=\left[\begin{matrix}1 & \Gamma_{j+1}z^{-1}\\ \Gamma^*_{j+1} & z^{-1}\end{matrix}\right]\left[\begin{matrix}G_{j}(z)\\ G_{j}^R(z)\end{matrix}\right]\tag{S.11}$

在这里插入图片描述

Recall that our goal is to derive a recursion that will take a sequence of autocorrelation values and generate the corresponding sequence of reflection coefficients:
$\left\{r_{x}(0), r_{x}(1), \ldots, r_{x}(p)\right\} \stackrel{\text {Schur}}{\longrightarrow}\left\{\Gamma_{1}, \Gamma_{2}, \ldots, \Gamma_{p}, \epsilon_{p}\right\}$
Therefore, we need to derive a method to find the reflection coefficient $\Gamma_{j+1}$ from $g_j(k)$ and $g_j^R(k)$ . Since $g_{j+1}(j+1)=0$ , evaluating $(S . 10)$ for $k = j + 1$ we have
$\Gamma_{j+1}=-\frac{g_j(j+1)}{g_j^R(j)}.$

在这里插入图片描述

Additional insight into the operation of the Schur recursion may be gained if it is formulated using vector notation.

From $(S . 9)$ :
$\left[\begin{matrix}\mathbf g_j^T\\(\mathbf g_j^R)^T\end{matrix}\right]=\left[\begin{matrix}\epsilon_j & \cdots & 0 & 0& g_j(j+1)& \cdots& g_j(p)\\0 & \cdots & 0 & g_j^R(j)& g_j^R(j+1)& \cdots& g^R_j(p)\end{matrix}\right]\tag{S.9}$

Shift the second row of the matrix to the right by one with $g_j^R(-1)$ entering as the first element of the second row (this corresponds to a delay of $g_j^R(k)$ by one or a multiplication of $G_j^R (z)$ by $z^{- 1}$ ).
$\left[\begin{matrix}\epsilon_j & \cdots & 0 & 0& g_j(j+1)& \cdots& g_j(p)\\g_j^R(-1) & \cdots & 0 &0 & g_j^R(j)& \cdots& g^R_j(p-1)\end{matrix}\right]$
(Note that $g_j^R(-1)=g_j(j+1)=\gamma_j$ and $g_j^R(j)=\epsilon_j$ ). After this shift, the ratio of the two terms in column number ( $j + 2$ ) is equal to $-\Gamma_{j+1}$ .
Applying the relationship in $(S . 10)$ , multiplying the shifted matrix by $\boldsymbol{\theta}_{j+1}=\left[\begin{matrix}1 & \Gamma_{j+1}\\ \Gamma_{j+1}^* & 1\end{matrix}\right]$ :
$\begin{aligned} \left[\begin{matrix}\mathbf g_{j+1}^T\\(\mathbf g_{j+1}^R)^T\end{matrix}\right]&=\left[\begin{matrix}1 & \Gamma_{j+1}\\ \Gamma_{j+1}^* & 1\end{matrix}\right]\left[\begin{matrix}\epsilon_j & \cdots & 0& g_j(j+1)& \cdots& g_j(p)\\g_j^R(-1) & \cdots &0 & g_j^R(j)& \cdots& g^R_j(p-1)\end{matrix}\right]\\ &=\left[\begin{matrix}\epsilon_{j+1} & \cdots & 0 & 0& g_{j+1}(j+2)&\cdots& g_{j+1}(p)\\0 & \cdots & 0 & g_{j+1}^R(j+1)& g^R_{j+1}(j+2)&\cdots& g^R_j(p-1)\end{matrix}\right] \end{aligned}$

This completes one step of the recursion.

Note that the first column of the matrix $(S . 9)$ never enters into the calculations, we may suppress the evaluation of these entries by initializing the first element in the first column to zero and by bringing in a zero into the second row each time that is shifted to the right. Since $g_j^R(j)=\epsilon_j, g_j^R(-1)=g_j(j+1)$ , no information is discarded with this simplification.

In summary, the steps described above lead to the following matrix formulation of the Schur recursion. Beginning with
$\mathbf{G}_{0}=\left[\begin{array}{ccccc} 0 & r_{x}(1) & r_{x}(2) & \cdots & r_{x}(p) \\ r_{x}(0) & r_{x}(1) & r_{x}(2) & \cdots & r_{x}(p) \end{array}\right]$
which is referred to as the generator matrix, a new matrix, $\widetilde{\mathbf{G}}_{0},$ is formed by shifting the second row of $\mathbf{G}_{0}$ to the right by one
$\widetilde{\mathbf{G}}_{0}=\left[\begin{array}{lllll} 0 & r_{x}(1) & r_{x}(2) & \cdots & r_{x}(p) \\ 0 & r_{x}(0) & r_{x}(1) & \cdots & r_{x}(p-1) \end{array}\right]$
Setting $\Gamma_{1}$ equal to the negative of the ratio of the two terms in the second column of $\widetilde{\mathbf{G}}_{0}$ , we then form the matrix
$\Theta_{1}=\left[\begin{array}{cc} 1 & \Gamma_{1} \\ \Gamma_{1}^{*} & 1 \end{array}\right]$
and evaluate $\mathbf{G}_{1}$ as follows
$\mathbf{G}_{1}=\boldsymbol{\Theta}_{1} \tilde{\mathbf{G}}_{0}=\left[\begin{array}{ccccc} 0 & 0 & g_{1}(2) & \cdots & g_{1}(p) \\ 0 & g_{1}^{R}(1) & g_{1}^{R}(2) & \cdots & g_{1}^{R}(p) \end{array}\right]$
The recursion then repeats these three steps where, in general, at the $j$ th step we

Shift the second row of $\mathbf{G}_{j}$ to the right by one,
Multiply $\widetilde{\mathbf{G}}_{j}$ by $\Theta_{j+1}$ to form $\mathbf{G}_{j+1}$ .

在这里插入图片描述

( $H$ in slides represents $\widetilde G$ in book)

Cholesky Factorization of a Toeplitz Matrix

Definition: For a positive matrix $C$ , the Cholesky factorization is a factorization as
$C=LDL^T,\quad L:\text{lower triangular, }D:\text{diagonal}$
The Levinson recursion or Schur algorithm provides a factorization of the Toeplitz matrix $R_x$ as follows:

Define the matrix $A_{p}$ in terms of the solutions of the Yule-Walker equations of order $0$ until $p :$
$\mathbf{A}_{p}=\left[\begin{array}{c|c|c|c|c} 1 & a_{1}(1) & a_{2}(2) & \vdots & a_{p}(p) \\ ~ & 1 & a_{2}(1) & \vdots & \vdots\\ ~ & ~ & 1 & \vdots & a_{p}(2) \\ ~ & ~ & ~ & 1& a_{p}(1) \\ ~ & ~ & ~ &~ & 1 \end{array}\right]$
Recall that for any order $p$ the Yule-Walker equation on the reverse sequence is
$\left[\begin{array}{ccccc} r_{x}(0) & r_{x}(1) & r_{x}(2) & \cdots & r_{x}(p) \\ r_{x}(1) & r_{x}(0) & r_{x}(1) & \cdots & r_{x}(p-1) \\ r_{x}(2) & r_{x}(1) & r_{x}(0) & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & r_{x}(1) \\ r_{x}(p) & r_{x}(p-1) & \cdots & r_{x}(1) & r_{x}(0) \end{array}\right]\left[\begin{array}{c} a_{p}(p) \\ \vdots \\ a_{p}(2) \\ a_{p}(1) \\ 1 \end{array}\right]=\left[\begin{array}{c} 0 \\ 0 \\ \vdots \\ 0 \\ \epsilon_{p} \end{array}\right]$
It follows that
$\left[\begin{array}{ccccc} r_{x}(0) & r_{x}(1) & r_{x}(2) & \cdots & r_{x}(p) \\ r_{x}(1) & r_{x}(0) & r_{x}(1) & \cdots & r_{x}(p-1) \\ r_{x}(2) & r_{x}(1) & r_{x}(0) & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & r_{x}(1) \\ r_{x}(p) & r_{x}(p-1) & \cdots & r_{x}(1) & r_{x}(0) \end{array}\right]\left[\begin{array}{c|c|c|c|c} 1 & a_{1}(1) & a_{2}(2) & \vdots & a_{p}(p) \\ ~ & 1 & a_{2}(1) & \vdots & \vdots\\ ~ & ~ & 1 & \vdots & a_{p}(2) \\ ~ & ~ & ~ & 1& a_{p}(1) \\ ~ & ~ & ~ &~ & 1 \end{array}\right]=\left[\begin{array}{c|c|c|c|c} \epsilon_0 & & & & \\ *& \epsilon_1 & & & \\ *& * & \epsilon_2 & & \\ *& * & * & *& \\ \vdots & \vdots & \vdots &\vdots & \epsilon_p \end{array}\right]\\ \Leftrightarrow \quad R_{x} A_{p}=E$
Consider now ${A}_{p}^{T} {R}_{x} {A}_{p} .$ Because ${A}_{p}^{T}$ and ${R}_{x} {A}_{p}$ are lower triangular matrices, it must be lower. But it is also a symmetric matrix. Hence it is a diagonal matrix. The entries on the main diagonal are seen to be $\left\{\epsilon_{0}, \cdots, \epsilon_{p}\right\}$ .
We thus found
${A}_{p}^{T} {R}_{x} {A}_{p}=D_{p} \quad \text { (diagonal) } \quad \Rightarrow \quad {R}_{x}={A}_{p}^{-T} {D}_{p} {A}_{p}^{-1}, \quad {R}_{x}^{-1}={A}_{p} {D}_{p}^{-1} {A}_{p}^{T}$
These are Cholesky factorizations of ${R}_{x}$ and ${R}_{x}^{-1}$

The Step-Up and Step-Down Recursions

Step-Up

$\left\{\Gamma_{1}, \Gamma_{2}, \ldots, \Gamma_{p}, \epsilon_{p}\right\} \stackrel{S t e p-u p}{\longrightarrow}\left\{a_{p}(1), a_{p}(2), \ldots, a_{p}(p), b(0)\right\}$

The Levinson order-update equation given in $(L F . 3)$ is a recursion for deriving the filter coefficients $a_{p}(k)$ from the reflection coefficients, $\Gamma_{i}$ . Specifically, since
$a_{j+1}(i)=a_{j}(i)+\Gamma_{j+1} a_{j}^{*}(j-i+1)$
then the filter coefficients $a_{j+1}(i)$ may be easily found from $a_{j}(i)$ and $\Gamma_{j+1} .$ The recursion is initialized by setting $a_{0}(0)=1$ and, after the coefficients $a_{p}(k)$ have been determined, the recursion is completed by setting $b(0)=\sqrt{\epsilon_{p}} .$

In matrix form, we have
$\left[\begin{array}{c} 1 \\ a_{j+1}(1) \\ a_{j+1}(2) \\ \vdots \\ a_{j+1}(j) \\ a_{j+1}(j+1) \end{array}\right]=\left[\begin{array}{c} 1 \\ a_{j}(1) \\ a_{j}(2) \\ \vdots \\ a_{j}(j) \\ 0 \end{array}\right]+\Gamma_{j+1}\left[\begin{array}{c} 0 \\ a_{j}^{*}(j) \\ a_{j}^{*}(j-1) \\ \vdots \\ a_{j}^{*}(1) \\ 1 \end{array}\right]$

在这里插入图片描述

Step-Down

$\left\{a_{p}(1), a_{p}(2), \ldots, a_{p}(p), b(0)\right\} \stackrel{\text {Step-down}}{\longrightarrow}\left\{\Gamma_{1}, \Gamma_{2}, \ldots, \Gamma_{p}, \epsilon_{p}\right\}$

Based on the fact that since
$\Gamma_j=a_j(j)$
the the reflection coefficients may be compute by running the Levinson-Durbin recursion backwards. Specifically, beginning with $\mathbf a_P$ we set $\Gamma_p=a_p(p)$ . Then, we recursively find each of the lower-order models, $\mathbf a_j$ , for $j=p-1,p-2,\cdots,1$ and set $\Gamma_j=a_j(j)$ as illustrated below:

From $(L F . 6)$ we have
$\left[\begin{array}{c} A_{j+1}(z) \\ A_{j+1}^{R}(z) \end{array}\right]=\left[\begin{array}{cc} 1 & \Gamma_{j+1} \\ \Gamma_{j+1}^* & 1 \end{array}\right]\left[\begin{array}{cc} 1 & 0 \\ 0 & z^{-1} \end{array}\right]\left[\begin{array}{c} A_{j}(z) \\ A_{j}^{R}(z) \end{array}\right]$
Hence
$\left[\begin{array}{c} A_{j}(z) \\ A_{j}^{R}(z) \end{array}\right]=\left[\begin{array}{cc} 1 & 0 \\ 0 & z \end{array}\right]\frac{1}{1-|\Gamma_{j+1}|^2} \left[\begin{array}{cc} 1 & -\Gamma_{j+1} \\ -\Gamma_{j+1}^* & 1 \end{array}\right] \left[\begin{array}{c} A_{j+1}(z) \\ A_{j+1}^{R}(z) \end{array}\right]$

在这里插入图片描述

Therefore
$A_{j}(z)=\frac{1}{1-\left|\Gamma_{j+1}\right|^{2}}\left[A_{j+1}(z)-\Gamma_{j+1} A_{j+1}^{R}(z)\right]$
Or by taking the reverse $z$ -transform,
$a_{j}(i)=\frac{1}{1-\left|\Gamma_{j+1}\right|^{2}}\left[a_{j+1}(i)-\Gamma_{j+1} a_{j+1}^{*}(j-i+1)\right]$
which is the step-down recursion. This recursion may also be written in vector form as follows
$\left[\begin{array}{c} a_{j}(1) \\ a_{j}(2) \\ \vdots \\ a_{j}(j) \end{array}\right]=\frac{1}{1-\left|\Gamma_{j+1}\right|^{2}}\left\{\left[\begin{array}{c} a_{j+1}(1) \\ a_{j+1}(2) \\ \vdots \\ a_{j+1}(j) \end{array}\right]-\Gamma_{j+1}\left[\begin{array}{c} a_{j+1}^{*}(j) \\ a_{j+1}^{*}(j-1) \\ \vdots \\ a_{j+1}^{*}(1) \end{array}\right]\right\}$
(Note that $\Gamma_{j+1}=a_{j+1}(j+1)$ and $\Gamma_j=a_j(j)$ )

在这里插入图片描述

Application: The Schur-Cohn Stability Test

This test is based on Property 2 (p.226), which states that the roots of a polynomial will lie inside the unit circle if and only if the magnitudes of the reflection coefficients are less than one.

Therefore, given a causal, linear shift-invariant filter with a rational system function,
$H(z)=\frac{B(z)}{A(z)}$
the filter may be tested for stability as follows:

First, the step-down recursion is applied to the coefficients of the denominator polynomial $A (z)$ to generate a reflection coefficient sequence, $\Gamma_j$ . The filter will then be stable if and only if all of the reflection coefficients are less than one in magnitude.