The Levinson-Durbin Recursion

Massachusetts Institute of Technology Department of Electrical Engineering and Computer Science 6.341: Discrete-Time Signal Processing OpenCourseWare 2006 Lecture 13 The Levinson-Durbin Recursion In the previous lecture we looked at all-pole signal modeling, linear prediction, and the stochastic inverse-whitening problem. Each scenario was related in concept to the problem of processing a signal s[n] by: � ⎤ � s[n] −� A1 1 − pk=1 ak z −k −� g[n], such that e[n] � g[n] − φ[n] � was minimized in some sense. In the deterministic case, ⎤we chose to minimize � = n e2 [n], � and in the stochastic case, we chose to minimize � = E e2 [n] . The causal solution for the parameters ak was found in each case to be of the form � ⎡ ⎡� ⎡ � �s [0] �s [1] a1 �s [1] · · · �s [p − 1] ⎥ �s [1] � ⎥ � ⎥ · · · �s [p − 2]� �s [0] ⎥ � ⎥a2 � ⎥�s [2]� = (1) ⎥ � ⎥ �, � ⎥ . . .. .. . .. .. ⎦ � ⎦ .. � ⎦ .. � . . . �s [p − 1] �s [p − 2] · · · �s [0] ap �s [p] or equivalently, Tp εp = r p . (2) In this lecture, we address the issue of how exactly to solve for the parameters a k (or equivalently for the vector εp ). The fact that εp can be found simply by calculating ε p = Tp−1 rp is worth a mention, but the focus of this lecture is on an elegant technique for finding ε p which is generally less computationally expensive than taking a matrix inversion and which is also recursive in the filter order p. This method for finding ε p is called the Levinson-Durbin recursion. In addition to being recursive in p, we’ll see that certain intermediate results from the recursion give coefficients for a lattice implementation of an LTI filter h[n] which implements our signal model as A �p . H(z) = S � (z) = 1 − k=1 ak z −k We’ll also see that the Levinson-Durbin recursion prescribes a similar implementation for the inverse filter 1/H(z). 1 Let’s first prove the Levinson-Durbin recursion. The basic idea of the recursion is to find the solution εp+1 for the (p + 1)st order case from the solution ε p for the pth order case. We’ll see that this hinges on the fact that Tp is a p × p Toeplitz matrix, that is, that T p is symmetric, and that all entries along a given diagonal are equal. This property of T p is a result of its definition: (Tp )ij � �s [i − j]. To begin the proof, consider what happens in the p = 1 case. Equation 2 becomes for p = 1 T1 ε1 = r 1 , or �s [0]a1 = �s [1], giving the solution a1 = �s [1] . �s [0] (3) The recursion will now be developed by evaluating Equation 2 for order p + 1: Tp+1 εp+1 = rp+1 . Because we will be dealing with a solution vector ε � for multiple orders α, we will adopt the (�) (�) notation ak to refer to the parameter ak for the αth-order model. (Note that a k is generally (m) not equal to ak .) Equation 2 is therefore represented in matrix form as � �s [0] �s [1] .. . �s [1] �s [0] .. . ··· ··· .. . ⎥ ⎥ ⎥ ⎥ ⎥ ⎦�s [p − 1] �s [p − 2] · · · �s [p] �s [p − 1] · · · ⎡ ⎡ � (p+1) ⎡ � a1 �s [1] �s [p − 1] �s [p] ⎥ (p+1) � ⎥ a2 �s [p − 2] �s [p − 1]� � ⎥ �s [2] � � �⎥ � � �⎥ . .. .. .. ⎥ .. � = ⎥ �. ⎥ �⎥ . . . � ⎥ � �⎥ � ⎦ (p+1) � �s [p] � �s [0] �s [1] � ⎦ap (p+1) �s [p + 1] �s [1] �s [0] a p+1 The matrix Tp+1 relates to Tp as � Tp+1 ⎥ ⎥ ⎥ =⎥ ⎥ ⎦ �s [p] �s [p − 1] .. . Tp �s [p] �s [p − 1] · · · �s [1] and the vector rp+1 relates to rp as rp+1 = � rp �s [p + 1] 2 � . �s [1] �s [0] ⎡ � � � �, � � Defining a new vector �p as rp upside-down, or � ⎡ �s [p] ⎥�s [p − 1]� ⎥ � �p = ⎥ �, ... ⎦ � �s [1] gives a more compact representation of T p+1 in terms of Tp : � � �p Tp . Tp+1 = (�p )T �s [0] Let’s now represent the (p + 1)st-order parameter vector ε p+1 in terms of the pth order vector εp , a correction term kp+1 and a correction vector �p as εp+1 = � εp 0 � + � �p kp+1 � . (4) Equation 2 for order p + 1 is therefore represented in terms of the pth order equation as � � �� � � �� � � �p Tp rp �p εp = + , 0 kp+1 �s [p + 1] (�p )T �s [0] which, sorting the equations out, implies Tp εp + Tp �p + �p kp+1 = rp (5) (�p )T εp + (�p )T �p + �s [0]kp+1 = �s [p + 1]. (6) and Substituting Equation 2 for order p = 1 into Equation 5 then gives −1 = −�p . Tp �p kp+1 (7) Note now that since Tp is Toeplitz, the matrix realization of the causal Yule-Walker equations for order p (Equation 1) implies also ⎡ � ⎡ � (p) ⎡ � ap �s [p] �s [0] �s [1] · · · �s [p − 1] ⎥ (p) � ⎥ � ⎥ �s [1] �s [0] · · · �s [p − 2]� ⎥�s [p − 1]� � ⎥a ⎥ p−1 � ⎥ � = �, ⎥ � ⎥ . .. .. . .. . � .. .. � �⎥ ⎦ . . . ⎦ .. � ⎦ (p) �s [1] �s [p − 1] �s [p − 2] · · · �s [0] a 1 or Tp ρp = � p , 3 where ρp is εp upside-down. This then implies, after substituting in Equation 5, �p = −kp+1 ρp . (8) We’ll now use this and Equation 6 to find k p+1 and thus have all the pieces needed to complete the recursion. Pre-multiplying both sides of Equation 8 by (� p )T gives the scalar relation (�p )T �p = − (�p )T ρp kp+1 . Note, however, that the scalar relation (� p )T ρp = (rp )T εp also holds, since rp is a re-ordered version of �p in the same way that εp is a re-ordered version of ρp . This allows us to re-write the above equation as (�p )T �p = − (rp )T εp kp+1 , which implies from Equation 6 (�p )T εp − (rp )T εp kp+1 + �s [0]kp+1 = �s [p + 1], or kp+1 = �s [p + 1] − (�p )T εp . �s [0] − (rp )T εp (9) Equations 4, 8, and 9 therefore define a recursion formula for finding the (p + 1)st-order model (p+1) parameters ak in terms of the previously-obtained pth-order solution. Equation 3 gives a (1) closed-form solution for a1 , which provides a starting point for the recursion. Since this ends a considerable algebraic detour, we’ll now summarize what was concluded. 4 The Levinson-Durbin recursion • Problem statement: for order p, solve • Definitions: ⎣ ⎢ (p) (p) (p) T εp = a1 , a 2 , . . . , a p ⎢ ⎣ (p) (p) (p) T ρp = ap , ap−1 , . . . , a1 (p) k=1 ak �s [i − k] = �s [i] ⎡ ap (p) ap−1 .. . �p rp = [�s [1], �s [2], . . . , �s [p]]T �p = [�s [p], �s [p − 1], . . . , �s [1]]T (1) • Solution for p = 1: a1 = �s [1] �s [0] • Recursion: kp+1 = �s [p + 1] − (�p )T εp �s [0] − (rp )T εp �p = −kp+1 ρp εp+1 = � εp 0 � + � �p kp+1 � � (p) a1 (p) a2 .. . ⎥ ⎥ ⎥ =⎥ ⎥ ⎥ (p) ⎦ ap 0 � (p) � ⎥ � ⎥ � ⎥ � ⎥ � − kp+1 ⎥ � ⎥ (p) � ⎦ a 1 −1 ⎡ � � � � � � � In parting, note that when these coefficients k � are used to implement an all-pole lattice filter with reflection coefficients k� , the lattice filter’s response is described by H(z) = S � (z) = A 1− �p m=1 am z −m . In other words, the impulse response of an all-pole lattice filter with reflection coefficients k � , α = 1, . . . , p, as prescribed by the Levinson-Durbin recursion, is the pth-order model of our original signal s[n]. (Likewise, the corresponding all-zero lattice filter implements 1/S � (z).) Note further that it is straightforward to determine the stability of an all-pole filter when im­ plemented in lattice form. Specifically, if the coefficients k � all have magnitude < 1, the all-pole filter is guaranteed to be stable. (It is left as an exercise to gain intuition for why this is so.) This stability criterion is of particular interest in applications where the all-pole filter is 5 implemented with finite-precision coefficients, since determining the stability of a comparable direct-form all-pole implementation requires the generally more-expensive process of determin­ ing all of the filter’s pole locations. In addition to the computational benefits gained from the Levinson-Durbin recursion, its connection with the lattice structure (and the structure’s asso­ ciated stability metric) enables a wide array of applications, including real-time speech coders and synthesizers. 6