1 / 44

Lecture 12

Lecture 12. Filter Theory. Review of last lecture output (“response”) of a linear system can be calculated by convolving its input (“forcing”) with its impulse response. Past history important. Output: temperature q (t).

trudy
Télécharger la présentation

Lecture 12

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Lecture 12 Filter Theory

  2. Review of last lectureoutput (“response”) of a linear systemcan be calculated byconvolving its input (“forcing”)with its impulse response

  3. Past history important Output: temperature q(t) “physics-based” response of temperature to an impulse in heating, g(t-t) Steel plate Input: heat h(t) convolution q(t) = -tg(t-t) h(t) dt

  4. h(t) h(t) t t t 0 q(t)=g(t) t amplitude h(t) h(t)g(t-t) q(t) t t 0 q(t) = -tg(t-t) h(t) dt

  5. Mathematical equivalent ways to write the convolution q(t) = -tg(t-t) h(t) dt or alternatively q(t) = 0g(t) h(t-t) dt h(t) is “forward in time” g(t) is “forward in time”

  6. h(t) hk t tk Dt approximation of functions as time-series hk=h[kDt] with k= … -2, -1, 0, 1, 2, … q(t) = -tg(t-t) h(t) dt  qk = Dt Sp=-k gk-p hp q(t) = 0g(t) h(t-t) dt  qk = Dt Sp=0 gp hk-p or alternatively

  7. q0 q1 … qN g0 0 0 0 0 0 g1 g0 0 0 0 0 … gN … g3 g2 g1 g0 h0 h1 … hN = Dt q0 q1 … qN h0 0 0 0 0 0 h1 h0 0 0 0 0 … hN … h3 h2 h1 h0 g0 g1 … gN = Dt Matrix formulations q = Gh and q = Gg

  8. End of review … q(t) = -tg(t-t) h(t) dt g(t) has a special meaning in terms of the impusle response of a dynamical system y(t) = -tf(t-t) x(t) dt A generic way of relating two functions x(t) and y(t)

  9. A generic way to construct a function of time y(t) is obtained from x(t) by convolving by filter f(t) y(t) = -tf(t-t) x(t) dt input filter output

  10. Similarly, a generic way to construct a time-series yk is obtained from xk by convolving by filter fk yk = Sp=-k fk-p xp input “digital” filter output

  11. By the way … The convolution is often written shorthand asy(t) = f(t)*x(t)y = f*x version using functions version using time-series Sure beats typing in integrals and summations …

  12. By the way …Here’s how to do convolution by hand x=[x0, x1, x2, x3, x4, …]T and y=[y0, y1, y2, y3, y4, …]T Reverse on time-series, line them up as shown, and multiply rows. This is first element of x*y x0, x1, x2, x3, x4, …  … y4, y3, y2, y1, y0 x0y0

  13. Slide to increase the overlap by one, multiply rows and add products. This is the second element x0, x1, x2, x3, x4, …   … y4, y3, y2, y1, y0 x0y1+x1y0 Slide again, multiply and add. This is the third element x0, x1, x2, x3, x4, …    … y4, y3, y2, y1, y0 x0y2+x1y1+x2y0 Repeat until time-series no longer overlap

  14. Filter Theory A collection of techniques used to design useful digital filters and to understand their behavior

  15. Sample filter theory question Given a filter fIs there a filter – call it finv – that undoesy = f*xso thatx = finv*y?

  16. A crazy but extremely important insight Convolution of time-series is identical to multiplication of polynomials

  17. computed via the convolution formula Convolution of time-series x=[x0, x1, x2, x3, x4, …]T f=[f0, f1, f2, f3, f4, …]T f*x=[f0x0, (f0x1+ f1x0), (f0x2+ f1x1 + f2x0), …]T Multiplication of polynomials x(z) = x0 + x1z + x2z2 + x3z3 + x4z4 + … f(z) = f0 + f1z + f2z2 + f3z3 + f4z4 + … f(z)x(z) = f0x0 + (f0x1+ f1x0)z + (f0x2+ f1x1 + f2x0)z2 … computed via polynomial multiplication

  18. Yes, the elements of f*x exactly match the coefficients of f(z)x(z) f*x = [f0x0, (f0x1+ f1x0), (f0x2+ f1x1 + f2x0), …]T f(z)x(z) = f0x0 + (f0x1+ f1x0)z + (f0x2+ f1x1 + f2x0)z2 …

  19. the z-transformturn a timeseries into a polynomialand vice versa time-series x=[x0, x1, x2, x3, x4, …]T polynomial x(z) = x0 + x1z + x2z2 + x3z3 + x4z4 + … Z-transform

  20. And why would we ever want to do such a thing?????????

  21. Because we know a lotabout polynomials!(or at least some people do …)

  22. Example:a filter or length Nf = [f0, f1, f2, f3, … fN-1]Tcorresponds to a polynomialof degree N-1 f(z) = f0 + f1z + f2z2 + f3z3 + f4z4 + …

  23. a polynomial of degree N-1can be factored into N-1 binomials each involving one of its roots, ri f(z) = (z-r1)(z-r2)…(z-rN-1) Note, by the way, that the roots are in general complex numbers, and that for real filters, they occur in conjugate pairs …

  24. But multiplying a sequence of binomials corresponds to convolving a sequence of length-2 filtersso any filter of length N can be written as a cascade of N-1 length-2 filtersf = [f0, f1, f2, f3, … fN-1]T = [-r1, 1]T* [-r2, 1]T*…* [-rN-1, 1]T

  25. Now what about the question … Given a filter fIs there a filter – call it finv – that undoesy = f*xso thatx = finv*yLet’s consider the simple casef = [1, -f1]T

  26. the filter f = [1, -f1]Tcorresponds to the polynomial f(z)=1-f1zThe operation y=f*xcorresponds to y(z)=f(z)x(z)so x(z)=y(z)/f(z)Thus the z-transform of the filter finv is the rational function 1/(1-f1z)

  27. But the rational function 1/(1-f1z)has Taylor series expansion1 + f1z + f12z2 + f13z3 + …Thusfinv = [1, f1, f12, f13, … ]T

  28. We learn a few important general things from this analysis

  29. f = [1, -f1]Tfinv = [1, f1, f12, f13, …]T#1The length of an inverse filter can be infinite, even though the length of the original filter is finite

  30. f = [1, -f1]Tfinv = [1, f1, f12, f13, …]T#2The inverse filter only exists when |f1|<1, for otherwise the elements of finv grow without bound

  31. since the general casef(z) = (z-r1)(z-r2)…(z-rN-1)can be rewrittenf(z) = (-r1)(1-r1-1z)  … (-rN-1)(1-rN-1-1z) #3In the general case, an inverse filter only exists when |ri-1|<1 or |ri|>1 That is, when all the roots of f(z) lie outside the “unit circle in the complex z-plane”(that just means |ri|>1)

  32. Lingoany filter that can be writtenas a cascade oflength-two filtersf = [f0, f1, f2, f3, … fN-1]T = c [1, -a1]T* [1, -a2]T*…* [1, -aN-1]T each of which satisfy |ai|<1is said to be minimum phase

  33. An interesting tidbit Multiplying x(z) by z is equivalent to shifting x by one sample to the right so f=[0, 1]T is the unit-delay filter x = [x0, x1, x2, x3, … xN-1]Tx(z) = x0 + x1z + x2z2 + x3z3 + x4z4 + …f(z)=zzx(z) = x0z + x1z2 + x2z3 + x3z4 + x4z5 + …f*x = [0, x0, x1, x2, x3, …]T

  34. Recursive filters Normal convolution is a lot of work when the filter length, N, is large yk = Sp=0N-1 fp xk-p N multiplications and additions … … at each of N time steps at which you want to know yk … so N2 overall

  35. A variant on the filtering idea is to use the y’s that you’ve already computed in the formula … yk = Sp=0N hp xk-p -Sp=1M dp yk-p this part “normal” filtering Here’s the recursion. It uses only up to yk-1 Thus past values of y are being used directly in the computation of the present value of y

  36. If we define d0=1, then the equation can be reqritten Sp=0M dp yk-p = Sp=0N hp xk-p or d*y = h*x Comparing to the usual y=f*x, we see that f = dinv*h

  37. In other words, might it be possible to choose two short filtersd and hso thatdinv*his a reasonably good approximation to amuch longer filter f?

  38. A Fairly Trivial Examplesuppose we wanted to match theinfinitely long filter f=[1, 1/2, 1/4,1/8,1/16, …]T this filter a smoothing filter in the sense that each y is an average of mostly the last few x’s

  39. Let’s tryd=[1, -d1]T and h=[1,0]T as we derived previouslydinv=[1, d1, d12,d13,d14, …]T and sodinvh=[1, d1, d12,d13,d14, …]Tthus the choice d1=1/2 matchesf=[1, 1/2, 1/4,1/8,1/16, …]T

  40. The recursive formula is thus yk = Sp=0N hp xk-p -Sp=1M dp yk-p yk = xk +½ yk-1 which involves very little computation

  41. Exemplary MatLab Code y=zeros(N,1); y(1)=x(1); for p = [2:N] y(p) = x(p) + 0.5*y(p-1); end

  42. x(t) time, t Red: normal filtering y(t) Black: recursive filtering time, t

  43. Note that the choice of d1 controls the amount of smoothing(the bigger the d1 the more the smoothing)so this filter is actually fairly flexible

  44. x(t) time, t y(t) d1=0.5 time, t y(t) d1=0.9 time, t

More Related