Hi,
I've a filter expressed as a parallel bank of (~50) second order sections.
Each second order section represents a high-Q resonant feature. I need to
invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} s.t. x
= F^{-1}y. (F is invertible.)
I'm stuck. Round-off precludes conversion to transfer function form; a
literature search (such as I'm able to do: I'm not trained in this area)
has not turned-up anything that has suggested a way forward. I'd appreciate
any advice from the experts on this list.
Thanks in advance,
Sam
Lee Samuel Finn
Professor of Physics; Astronomy & Astrophysics
104 Davey Laboratroy
The Pennsylvania State University
University Park, PA 16802
Ph: +1(814)863-9598
Fx: +1(814)863-9608
invert filter bank of parallel second order sections?
Started by ●July 11, 2010
Reply by ●July 11, 20102010-07-11
On 11 Jul, 23:30, "lsfinn" <lsfinn@n_o_s_p_a_m.psu.edu> wrote:> Hi, > > I've a filter expressed as a parallel bank of (~50) second order sections. > Each second order section represents a high-Q resonant feature. I need to > invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} s.t. x > = F^{-1}y. (F is invertible.) > > I'm stuck. Round-off precludes conversion to transfer function form; a > literature search (such as I'm able to do: I'm not trained in this area) > has not turned-up anything that has suggested a way forward. I'd appreciate > any advice from the experts on this list.Won't claim any expert status, but apart form that: First of all, what you want might not be possible. The deciding factor is whether the second-order sections (SOS) are minimum phase. Denote the transfer function of one section as H(z) = B(z)/A(z). The expression for the inverted filter will then be H'(z) = A(z)/B(z). The usual condition for this filter to be causal and stable is that all the roots of the polynomial B(z) are inside the unit circle in Z domain. Inverting each SOS this way works for a cascade of SOS'es. I can't tell off the top of my head if it would work for a parallel bank of SOS'es; maybe some of the true experts will chime in on that. Rune
Reply by ●July 11, 20102010-07-11
lsfinn wrote:> Hi, > > I've a filter expressed as a parallel bank of (~50) second order sections. > Each second order section represents a high-Q resonant feature. I need to > invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} s.t. x > = F^{-1}y. (F is invertible.)You might want to ask RBJ; this is his area. My $0.02: The straightforward way would be collapse the filter into transfer function, factorize it and then invert. Using multiple precision math. There could be many ways to approximate the inverse or to solve for the inverse as numeric optimization problem. Vladimir Vassilevsky DSP and Mixed Signal Design Consultant http://www.abvolt.com
Reply by ●July 11, 20102010-07-11
On 07/11/2010 03:17 PM, Rune Allnor wrote:> On 11 Jul, 23:30, "lsfinn"<lsfinn@n_o_s_p_a_m.psu.edu> wrote: >> Hi, >> >> I've a filter expressed as a parallel bank of (~50) second order sections. >> Each second order section represents a high-Q resonant feature. I need to >> invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} s.t. x >> = F^{-1}y. (F is invertible.) >> >> I'm stuck. Round-off precludes conversion to transfer function form; a >> literature search (such as I'm able to do: I'm not trained in this area) >> has not turned-up anything that has suggested a way forward. I'd appreciate >> any advice from the experts on this list. > > Won't claim any expert status, but apart form that: > > First of all, what you want might not be possible. The deciding > factor is whether the second-order sections (SOS) are minimum > phase. Denote the transfer function of one section as > > H(z) = B(z)/A(z). > > The expression for the inverted filter will then be > > H'(z) = A(z)/B(z). > > The usual condition for this filter to be causal and stable > is that all the roots of the polynomial B(z) are inside the > unit circle in Z domain. > > Inverting each SOS this way works for a cascade of SOS'es. > I can't tell off the top of my head if it would work for > a parallel bank of SOS'es; maybe some of the true experts > will chime in on that.Alas, this is one of those "I can disprove that hypothesis in one step" things. Consider the filter H(z) = (z - a) / (z - b) + c -- it's parallel, and each leg is minimum phase (or close enough), right? H(z) = ((1 + c)*z - a - b*c) / (z - b). Just find a, b, and c such that (a + b*c) / (1 + c) > 1 Now you've found a minimum-phase system from your parallel filter bank. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Do you need to implement control loops in software? "Applied Control Theory for Embedded Systems" was written for you. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●July 11, 20102010-07-11
On 07/11/2010 03:34 PM, Vladimir Vassilevsky wrote:> > > lsfinn wrote: > >> Hi, >> I've a filter expressed as a parallel bank of (~50) second order >> sections. >> Each second order section represents a high-Q resonant feature. I need to >> invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} >> s.t. x >> = F^{-1}y. (F is invertible.) > > You might want to ask RBJ; this is his area. > > My $0.02: > > The straightforward way would be collapse the filter into transfer > function, factorize it and then invert. Using multiple precision math.I second that.> There could be many ways to approximate the inverse or to solve for the > inverse as numeric optimization problem.For that matter, if what you're _really_ trying to do is to take the filter's output signal and deduce it's input signal, then you may well do better with a Kalman (or Wiener, IOW a "time-invariant Kalman"). This would let you take any noise into account (physical or numeric) that's added in _after_ the filtering. I.e., if you have the signal y = h(x) + n, where x is the input signal to your filter bank and n is a noise signal, and you want the best estimate of y, then an inverse filter bank may _not_ be the best filter at all. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Do you need to implement control loops in software? "Applied Control Theory for Embedded Systems" was written for you. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●July 12, 20102010-07-12
On Jul 11, 6:34�pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote:> lsfinn wrote: > > Hi, > > > I've a filter expressed as a parallel bank of (~50) second order sections. > > Each second order section represents a high-Q resonant feature. I need to > > invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} s.t. x > > = F^{-1}y. (F is invertible.) > > You might want to ask RBJ; this is his area.i don't think necessarily so.> > My $0.02: > > The straightforward way would be collapse the filter into transfer > function, factorize it and then invert. Using multiple precision math. > > There could be many ways to approximate the inverse or to solve for the > inverse as numeric optimization problem.Vlad, i can't add a penny to what you or Tim have said. if Sam has access to both x and y of his filter bank, he could try some sorta adaptive filter in series where the input is y and the desired signal is x. r b-j
Reply by ●July 12, 20102010-07-12
On 07/11/2010 08:35 PM, robert bristow-johnson wrote:> On Jul 11, 6:34 pm, Vladimir Vassilevsky<nos...@nowhere.com> wrote: >> lsfinn wrote: >>> Hi, >> >>> I've a filter expressed as a parallel bank of (~50) second order sections. >>> Each second order section represents a high-Q resonant feature. I need to >>> invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} s.t. x >>> = F^{-1}y. (F is invertible.) >> >> You might want to ask RBJ; this is his area. > > i don't think necessarily so. > >> >> My $0.02: >> >> The straightforward way would be collapse the filter into transfer >> function, factorize it and then invert. Using multiple precision math. >> >> There could be many ways to approximate the inverse or to solve for the >> inverse as numeric optimization problem. > > Vlad, i can't add a penny to what you or Tim have said. > > if Sam has access to both x and y of his filter bank, he could try > some sorta adaptive filter in series where the input is y and the > desired signal is x.Yes, I think we need to know more. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Do you need to implement control loops in software? "Applied Control Theory for Embedded Systems" was written for you. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●July 12, 20102010-07-12
On Sun, 11 Jul 2010 15:17:32 -0700 (PDT), Rune Allnor <allnor@tele.ntnu.no> wrote:>On 11 Jul, 23:30, "lsfinn" <lsfinn@n_o_s_p_a_m.psu.edu> wrote: >> Hi, >> >> I've a filter expressed as a parallel bank of (~50) second order sections. >> Each second order section represents a high-Q resonant feature. I need to >> invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} s.t. x >> = F^{-1}y. (F is invertible.)...>Denote the transfer function of one section as > >H(z) = B(z)/A(z). > >The expression for the inverted filter will then be > >H'(z) = A(z)/B(z).The parallel bank of second order sections is the partial fraction decomposition of the full transfer function. The transfer function therefore has the same poles as the individual filters. But it does not have the same zeroes. Rudimentary example: X Xz^-1 --- = ------- (X,x arbitrary) has a zero at infinity. z+x 1+xz^-1 A B A(z+b) + B(z+a) (A+B)z + Ab + aB z + (Ab+aB)/(A+B) --- + --- = --------------- = ---------------- = (A+B)----------------- z+a z+b (z+a)(z+b) (z+a)(z+b) (z+a)(z+b) has a zero at -(Ab+aB)/(A+B). As the OP mentioned, numerical precision problems make computing the full transfer function impossible. However, knowing the poles of the full transfer function (which are the zeroes of the inverted transfer function) should give enough a priori information for a decent system identification. Think of a formulation of Prony's Method in which the poles are already known exactly http://cnx.org/content/m12762/latest/, equation 3. Greg
Reply by ●July 12, 20102010-07-12
robert bristow-johnson wrote:> On Jul 11, 6:34 pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote: > >>lsfinn wrote: >> >>>Hi, >> >>>I've a filter expressed as a parallel bank of (~50) second order sections. >>>Each second order section represents a high-Q resonant feature. I need to >>>invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} s.t. x >>>= F^{-1}y. (F is invertible.) >> >>You might want to ask RBJ; this is his area. > > i don't think necessarily so.I far as I remember you worked on the problem of decomposing a transfer function into the *parallel* set of biquads. It is not obvious to me what to do with numerator. Vladimir Vassilevsky DSP and Mixed Signal Design Consultant http://www.abvolt.com
Reply by ●July 12, 20102010-07-12
On Jul 12, 1:31�pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote:> robert bristow-johnson wrote: > > On Jul 11, 6:34 pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote: > > >>lsfinn wrote: > > >>>Hi, > > >>>I've a filter expressed as a parallel bank of (~50) second order sections. > >>>Each second order section represents a high-Q resonant feature. I need to > >>>invert the filter: i.e., if F = A1+A2+... and y = Fx, I need F^{-1} s.t. x > >>>= F^{-1}y. (F is invertible.) > > >>You might want to ask RBJ; this is his area. > > > i don't think necessarily so. > > I far as I remember you worked on the problem of decomposing a transfer > function into the *parallel* set of biquads.which is the opposite direction of what Sam is doing. also, i didn't get very far with that idea, but i should return to it.> It is not obvious to me what to do with numerator.well, if the original composite filter was known to be minimum phase, you switch numerator and denominator. then break it out again into partial fractions, if that is what Sam wants. r b-j
