10: Digital Filter

⊲ Structures

Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage

10: Digital Filter Structures

Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example Lattice Example Numerator Summary MATLAB routines

DSP and Digital Filters (2015-7227)

Structures: 10 – 1 / 19

Direct Forms 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads

Filter: H(z) = B(z) A(z) with input x[n] and output y[n] PM PN y[n] = k=0 b[k]x[n − k] − k=1 a[k]y[n − k]

Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage

Direct Form 1: • Direct implementation of difference equation 1 • Can view as B(z) followed by A(z)

⊲

Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

Direct forms use coefficients a[k] and b[k] directly

Direct Form II: 1 followed by B(z) • Implements A(z) • Saves on delays (= storage)

Lattice Example Numerator Summary MATLAB routines

DSP and Digital Filters (2015-7227)

Structures: 10 – 2 / 19

Transposition 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads

⊲

Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

Can • • • •

convert any block diagram into an equivalent transposed form: Reverse direction of each interconnection Reverse direction of each multiplier Change junctions to adders and vice-versa Interchange the input and output signals

Example: Direct form II → Direct Form IIt Would normally be drawn with input on the left Note: A valid block diagram must never have any feedback loops that don’t go through a delay (z −1 block).

Lattice Example Numerator Summary MATLAB routines

DSP and Digital Filters (2015-7227)

Structures: 10 – 3 / 19

State Space 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads

⊲

Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example Lattice Example Numerator Summary MATLAB routines

v[n] is a vector of delay element outputs Can write: v[n + 1] = Pv[n] + qx[n] y[n] = rT v[n] + sx[n]  T P, q, r , s is the state-space representation of the filter structure. The transfer function is given by: H(z) =

B(z) A(z)

=

det(zI−P+qrT ) det(zI−P) T

+s−1

The transposed form has P → P and q ↔ r Example: Direct Form IIt   −a[1] 1 P= −a[2] 0
T r = 1 0 From which H(z) =

DSP and Digital Filters (2015-7227)

q=



same H(z)

⇒

b[1] − b[0]a[1] b[2] − b[0]a[2]



s = b[0] b[0]z 2 +b[1]z+b[2] z 2 +a[1]z+a[2] Structures: 10 – 4 / 19

[State-Space → Transfer Function] [This is not examinable] We start by proving a useful formula which shows how the determinant of a matrix, A, changes when you add a rank-1 matrix, qrT , onto it. The formula is known as the Matrix Determinant Lemma. For any matrix A and column vectors q and r, we can write       T T −1 T T T 1 r 1+r A q 0 1 r 1 0 = . −1 0 A −q I −A q I 0 A + qrT It is easy to verify this by multiplying out the matrices. We now take the determinant of both sides making use of the result than a block triangular matrix is the product of the blocks along the diagonal (assuming they are all square). This gives the result we want:

det (A) × 1 + rT A−1 q = det A + qrT Now we take the z-transform of the state space equations v[n + 1] = Pv[n] + qx[n] y[n] = rT v[n] + sx[n]

−→

z−transform

zV = PV + qX Y = rT V + sX

The upper equation gives (zI − P) V = qX from which V = (zI − P)−1 qX and by substituting this det(zI−P+qrT ) −1 Y T + s − 1. in the lower equation, we get X = r (zI − P) q + s = det(zI−P) DSP and Digital Filters (2015-7227)

Structures: 10 – note 1 of slide 4

Precision Issues 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads

⊲

Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example Lattice Example Numerator Summary MATLAB routines

If all computations were exact, it would not make any difference which of the equivalent structures was used. However ... • Coefficient precision Coefficients are stored to finite precision and so are not exact. The filter actually implemented is therefore incorrect. • Arithmetic precision Arithmetic calculations are not exact. ◦ Worst case for arithmetic errors is when calculating the difference between two similar values: 1.23456789 − 1.23455678 = 0.00001111: 9 s.f. → 4 s.f. Arithmetic errors introduce noise that is then filtered by the transfer function between the point of noise creation and the output.

DSP and Digital Filters (2015-7227)

Structures: 10 – 5 / 19

Coefficient Sensitivity 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads

⊲

Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

The roots of high order polynomials can be very sensitive to small changes in coefficient values. Wilkinson’s polynomial: (famous example) Q20 f (x) = n=1 (x − n) = x20 − 210x19 + 20615x18 − . . . has roots well separated on the real axis.

Multiplying the coefficient of x19 by 1.000001 moves the roots a lot. “Speaking for myself I regard it as the most traumatic experience in my career as a numerical analyst”, James Wilkinson 1984

Lattice Example Numerator Summary MATLAB routines

5

5

0

0

-5

-5

0

5

10

15

20

25

0

5

10

15

20

25

Moral: Avoid using direct form for filters orders over about 10.

DSP and Digital Filters (2015-7227)

Structures: 10 – 6 / 19

Cascaded Biquads 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads

Avoid high order polynomials by factorizing into quadratic terms: B(z) A(z)

=

Q (1+bk,1 z−1 +bk,2 z−2 ) Q g (1+ak,1 z−1 +ak,2 z−2 ) =

⊲

Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

The term

g

1+bk,1 z −1 +bk,2 z −2 k=1 1+ak,1 z −1 +ak,2 z −2

QK

 M   N 
where K = max 2 , 2 . 1+bk,1 z −1 +bk,2 z −2 1+ak,1 z −1 +ak,2 z −2

is a biquad (bi-quadratic section).

We need to choose: (a) which poles to pair with which zeros in each biquad (b) how to order the biquads

Lattice Example Numerator Summary MATLAB routines

Direct Form II Transposed

DSP and Digital Filters (2015-7227)

Structures: 10 – 7 / 19

Pole-zero Pairing/Ordering

Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage

⊲

Example: Elliptic lowpass filter

1 0.5

2 pole pairs and 2 zero pairs need 2 biquads

0 -0.5 -1

Noise introduced in one biquad is amplified by all the subsequent ones:

-1

0 z

1

• Make the peak gain of each biquad as small as possible ◦ Pair poles with nearest zeros to get lowest peak gain begin with the pole nearest the unit circle ◦ Pairing with farthest zeros gives higher peak biquad gain

Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

• Poles near the unit circle have the highest peaks and introduce most noise so place them last in the chain

Lattice Example Numerator Summary MATLAB routines

Nearest

Gain (dB)

20 0 -20 -40 0

0.5

1

DSP and Digital Filters (2015-7227)

1.5 ω

2

2.5

3

Farthest

20 Gain (dB)

10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads

0 -20 -40 0

0.5

1

1.5 ω

2

2.5

3

Structures: 10 – 8 / 19

Linear Phase 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage

⊲

Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

Implementation can take advantage of any symmetry in the coefficients. Linear phase filters are always FIR and have symmetric (or, more rarely, antisymmetric) coefficients. PM H(z) = m=0 h[m]z −m h[M − m] = h[m]  M  − M P M2 −1
−m m−M 2 + m=0 h[m] z =h 2 z [m even] +z For M even, we only need M 2 + 1 multiplies instead of M + 1. We need M additions in each case.

Lattice Example Numerator Summary MATLAB routines

M = 6: For M odd (no central coefficient), we need

DSP and Digital Filters (2015-7227)

M +1 2

multiplies.

Structures: 10 – 9 / 19

Hardware Implementation 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage

⊲

Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

Software Implementation: All that matters is the total number of multiplies and adds Hardware Implementation: Delay elements (z −1 ) represent storage registers The maximum clock speed is limited by the number of sequential operations between registers Example: Symmetric Linear Phase Filter Direct form: Maximum sequential delay = 4a + m Transpose form: Maximum sequential delay = a + m , a and m are the delays of adder and multiplier respectively

Lattice Example Numerator Summary MATLAB routines

DSP and Digital Filters (2015-7227)

Structures: 10 – 10 / 19

Allpass Filters 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage

Allpass filters have mirror image numerator and denominator coefficients: b[n] = a[N − n] ⇔ jω ⇒ H(e ) ≡ 1∀ω

There are several efficient structures, e.g. • First Order: H(z) =

a[1]+z −1 1+a[1]z −1

⊲

Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

B(z) = z −N A(z −1 )

• Second Order: H(z) =

a[2]+a[1]z −1 +z −2 1+a[1]z −1 +a[2]z −2

Lattice Example Numerator Summary MATLAB routines

Allpass filters have a gain magnitude of 1 even with coefficient errors.

DSP and Digital Filters (2015-7227)

Structures: 10 – 11 / 19

Lattice Stage 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage

⊲

Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example Lattice Example Numerator Summary MATLAB routines

Suppose G is allpass: G(z) =

z −N A(z −1 ) A(z)

V (z) = X(z) − kGz −1 V (z) ⇒ V (z) =

1 1+kGz −1 X(z)

Y (z) = kV (z) + Gz Y (z) X(z)

=

−1

V (z) =

kA(z)+z −N −1 A(z −1 ) A(z)+kz −N −1 A(z −1 )

,

k+z −1 G 1+kGz −1 X(z)

z −(N +1) D(z −1 ) D(z)

Obtaining {d[n]} from {a[n]}:   n=0 1 d[n] = a[n] + ka[N + 1 − n] 1 ≤ n ≤ N   k n=N +1 Obtaining {a[n]} from {d[n]}: k = d[N + 1] a[n] = If G(z) is stable then

DSP and Digital Filters (2015-7227)

Y (z) X(z)

d[n]−kd[N +1−n] 1−k2

is stable if and only if |k| < 1 (see note) Structures: 10 – 12 / 19

[Proof of Stability Criterion] We want to show that if G(z) is a stable allpass filter then |k| < 1.

Y (z) X(z)

=

k+z −1 G(z) 1+kz −1 G(z)

is stable if and only if

We make use of a property of allpass filters (proved in a note in lecture 5) that if G(z) is a stable allpass filter, then |G(z)| T 1 according to whether |z| S 1. If z is a root of the denominator 1 + kz −1 G(z), then kz −1 G(z)

=

−1

⇒ |k| × |z −1 | × |G(z)|

=

1

⇒ |k|

=

|z| |G(z)|

It follows from the previously stated property of G(z) that |z| S 1 ⇔

DSP and Digital Filters (2015-7227)

|z| |G(z)|

S 1 ⇔ |k| S 1.

Structures: 10 – note 1 of slide 12

Example A(z) ↔ D(z) 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

⊲

Lattice Example Numerator Summary MATLAB routines

Suppose N = 3, k = 0.5 and A(z) = 1 + 4z −1 − 6z −2 + 10z −3

A(z) → D(z) z0 1

A(z) z −4 A(z −1 ) D(z) = A(z) + kz −4 A(z −1 )

1

z −1 4 10 9

z −2 −6 −6 −9

z −3 10 4 12

z −4 1 0.5

D(z) → A(z) D(z) k = d[N + 1] z −4 D(z −1 ) D(z) − kz −4 D(z −1 ) A(z) =

D(z)−kz −4 D(z −1 ) 1−k2

DSP and Digital Filters (2015-7227)

z0 1

z −1 9

z −2 −9

z −3 12

0.5 0.75

12 3

−9 −4.5

9 7.5

z −4 0.5 0.5 1 0

1

4

−6

10

0

Structures: 10 – 13 / 19

Allpass Lattice 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

⊲

Lattice Example Numerator Summary MATLAB routines

We can implement any allpass filter H(z) = with M stages:

z −M A(z −1 ) A(z)

as a lattice filter

• Initialize AM (z) = A(z) • Repeat for m = M : −1 : 1 ◦ k[m] = am [m] ◦ am−1 [n] =

am [n]−k[m]am [m−n] 1−k2 [m]

equivalently Am−1 (z) =

for 0 ≤ n ≤ m − 1

Am (z)−k[m]z −m Am (z −1 ) 1−k2 [m]

A(z) is stable iff |k[m]| < 1 for all m (good stability test)

DSP and Digital Filters (2015-7227)

Structures: 10 – 14 / 19

Lattice Filter

Label outputs um [n] and vm [n] and define Hm (z) =

Vm (z) Um (z)

=

z −m Am (z −1 ) Am (z)

From earlier slide: Um−1 (z) Um (z)

=

1 1+k[m]z −1 H

m−1 (z)

=

Am−1 (z) Am−1 (z)+k[m]z −m Am−1 (z −1 )

=

Am−1 (z) Am (z)

Hence: Um (z) X(z)

=

Am (z) A(z)

and

Vm (z) X(z)

=

Um (z) X(z)

×

Vm (z) Um (z)

=

z −m Am (z −1 ) A(z)

m (z) is of order m so you can create any numerator of order M by The numerator of VX(z) summing appropriate multiples of Vm (z): PM −m PM Am (z −1 ) m=0 cm z w[n] = m=0 cm vm [n] ⇒ W (z) = A(z)

DSP and Digital Filters (2015-7227)

Structures: 10 – 15 / 19

Lattice Example

A(z) = A3 (z) = 1 + 0.2z −1 − 0.23z −2 + 0.2z −3 • k[3] = 0.2 ⇒a2 [ ] =

[1, 0.2, −0.23]−0.2[0.2, −0.23, 0.2] 1−0.22

• k[2] = −0.281 ⇒a1 [ ] =

= [1, 0.256, −0.281]

[1, 0.256]+0.281[−0.281, 0.256] 1−0.2812

= [1, 0.357]

• k[1] = 0.357 ⇒a0 [ ] = 1 V0 (z) X(z)

=

1 −1 1+0.2z −0.23z −2 +0.2z −3

V1 (z) X(z)

=

0.357+z −1 1+0.2z −1 −0.23z −2 +0.2z −3

V2 (z) X(z)

=

−0.281+0.256z −1 +z −2 1+0.2z −1 −0.23z −2 +0.2z −3

V3 (z) X(z)

=

0.2−0.23z −1 +0.2z −2 +z −3 1+0.2z −1 −0.23z −2 +0.2z −3

Add together multiples of DSP and Digital Filters (2015-7227)

Vm (z) X(z)

to create an arbitrary

B(z) 1+0.2z −1 −0.23z −2 +0.2z −3 Structures: 10 – 16 / 19

Lattice Example Numerator PM

Form a new output signal as w[n] = m=0 cm vm [n] PM B(z) W (z) = m=0 cm Vm (z) = 1+0.2z−1 −0.23z −2 +0.2z −3 X(z) V0 (z) X(z)

=

1 −1 1+0.2z −0.23z −2 +0.2z −3

V1 (z) X(z)

=

0.357+z −1 1+0.2z −1 −0.23z −2 +0.2z −3

V3 (z) −0.281+0.256z +z 0.2−0.23z +0.2z +z = 1+0.2z −1 −0.23z −2 +0.2z −3 X(z) = 1+0.2z −1 −0.23z −2 +0.2z −3      c0 1 0.357 −0.281 0.2 b[0]  c1   b[1]   0 1 0.256 −0.23        = We have    c2  0 0 1 0.2 b[2] c3 0 0 0 1 b[3]    −1  c0 1 0.357 −0.281 0.2 b[0]  c1   0  b[1] 1 0.256 −0.23       Hence choose cm as  =    b[2] c2 0 0 1 0.2 c3 0 0 0 1 b[3] V2 (z) X(z)

DSP and Digital Filters (2015-7227)

−1

−2

−1

−2

−3

   

Structures: 10 – 17 / 19

Summary 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example Lattice Example Numerator Summary MATLAB routines

• Filter block diagrams ◦ Direct forms ◦ Transposition ◦ State space representation • Precision issues: coefficient error, arithmetic error ◦ cascaded biquads • Allpass filters ◦ first and second order sections • Lattice filters ◦ Arbitrary allpass response ◦ Arbitrary IIR response by summing intermediate outputs

⊲

DSP and Digital Filters (2015-7227)

For further details see Mitra: 8.

Structures: 10 – 18 / 19

MATLAB routines 10: Digital Filter Structures Direct Forms Transposition State Space Precision Issues Coefficient Sensitivity Cascaded Biquads Pole-zero Pairing/Ordering Linear Phase Hardware Implementation Allpass Filters Lattice Stage Example A(z) ↔ D(z) Allpass Lattice Lattice Filter Lattice Example

residuez tf2sos,sos2tf zp2sos,sos2zp zp2ss,ss2zp tf2ss,ss2tf poly

b(z −1) a(z −1 ) b(z −1) a(z −1 )

↔

→ k 1−prkkz−1 Q b0,l +b1,l z−1 +b2,l z−2

{zm , pk , g} ↔

P

l 1+a1,l z −1 +a2,l z −2 Q b0,l +b1,l z−1 +b2,l z−2 (l 1+a∈1,l z−1 +a2,l z−2 ′

{zm , pk , g} ↔ (

x = Ax + Bu y = Cx + Du

x′ = Ax + Bu ↔ y = Cx + Du poly(A) = det (zI−A)

b(z −1) a(z −1 )

Lattice Example Numerator Summary MATLAB routines

⊲

DSP and Digital Filters (2015-7227)

Structures: 10 – 19 / 19