440 likes | 589 Views
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).
E N D
Lecture 12 Filter Theory
Review of last lectureoutput (“response”) of a linear systemcan be calculated byconvolving its input (“forcing”)with its impulse response
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
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
Mathematical equivalent ways to write the convolution q(t) = -tg(t-t) h(t) dt or alternatively q(t) = 0g(t) h(t-t) dt h(t) is “forward in time” g(t) is “forward in time”
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) = 0g(t) h(t-t) dt qk = Dt Sp=0 gp hk-p or alternatively
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
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)
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
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
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 …
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
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
Filter Theory A collection of techniques used to design useful digital filters and to understand their behavior
Sample filter theory question Given a filter fIs there a filter – call it finv – that undoesy = f*xso thatx = finv*y?
A crazy but extremely important insight Convolution of time-series is identical to multiplication of polynomials
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
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 …
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
Because we know a lotabout polynomials!(or at least some people do …)
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 + …
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 …
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
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
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)
But the rational function 1/(1-f1z)has Taylor series expansion1 + f1z + f12z2 + f13z3 + …Thusfinv = [1, f1, f12, f13, … ]T
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
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
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)
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
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
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
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
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
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?
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
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
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
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
x(t) time, t Red: normal filtering y(t) Black: recursive filtering time, t
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
x(t) time, t y(t) d1=0.5 time, t y(t) d1=0.9 time, t