Hi, I have an 8th order allpass filter and want to implement it as a cascade of 2nd order ALLPASS sections. The usual 'split into SOS' method doesn't keep the allpass property, so I now have to manually look at the pole/zeros of the high order filter, pick the pairs that are 2nd order friendly, and rebuild the filter. Is there a better way to do this - maybe this factoring with allpass constraint has a 'name'. Thanks!
Factoring a high order allpass filter into cascade of 2nd order ALLPASS
Started by ●May 10, 2012
Reply by ●May 10, 20122012-05-10
On 5/10/12 12:23 PM, gretzteam wrote:> Hi, > I have an 8th order allpass filter and want to implement it as a cascade of > 2nd order ALLPASS sections. The usual 'split into SOS' method doesn't keep > the allpass property,turn the APF into LPF (ignoring the numerator), factor into Second Order Sections, then reflect each pole explicitly to the APF zero.> so I now have to manually look at the pole/zeros of > the high order filter, pick the pairs that are 2nd order friendly, and > rebuild the filter.do the same (rebuild the filter), but if you are guaranteed that it's APF to start with, knowledge of the zeros doesn't do anything for you. just factor the poles and reflect each pole to its corresponding zero, rather than trying to find the zero that matches the pole.> > Is there a better way to do this - maybe this factoring with allpass > constraint has a 'name'. > > Thanks!FWIW -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●May 10, 20122012-05-10
> >do the same (rebuild the filter), but if you are guaranteed that it's >APF to start with, knowledge of the zeros doesn't do anything for you. >just factor the poles and reflect each pole to its corresponding zero, >rather than trying to find the zero that matches the pole. >>> Thanks! > >FWIW >That works great, thanks! Definitely easier to code for different filter orders!
Reply by ●May 10, 20122012-05-10
On Thu, 10 May 2012 12:40:12 -0400, robert bristow-johnson wrote:> On 5/10/12 12:23 PM, gretzteam wrote: >> Hi, >> I have an 8th order allpass filter and want to implement it as a >> cascade of 2nd order ALLPASS sections. The usual 'split into SOS' >> method doesn't keep the allpass property, > > turn the APF into LPF (ignoring the numerator), factor into Second Order > Sections, then reflect each pole explicitly to the APF zero. > >> so I now have to manually look at the pole/zeros of the high order >> filter, pick the pairs that are 2nd order friendly, and rebuild the >> filter. > > do the same (rebuild the filter), but if you are guaranteed that it's > APF to start with, knowledge of the zeros doesn't do anything for you. > just factor the poles and reflect each pole to its corresponding zero, > rather than trying to find the zero that matches the pole. > > > >> Is there a better way to do this - maybe this factoring with allpass >> constraint has a 'name'.If there's anything that I'd add to Robert's advise, it is this: Your problem most likely stems from trying to factor high-degree polynomials. Polynomial roots are notoriously sensitive to coefficient variations, and get ever more so as the order gets higher. Yet IIR filters are sensitive to the _roots_ of their numerator and denominator -- not the actual coefficients of the overall transfer function. Which means that your problem with your gazzilionth-order numerator and gazzilionth-order denominator isn't that you are getting your desired pole positions but not your desired zeros -- it means that _both_ your poles and your zeros are screwed up, but they just happen to be screwed up differently. Using Robert's method is essentially saying "OK, I'm going to live with the wrong pole positions (and hence a phase response that isn't to my original design) as long as my amplitude stays flat. So _my_ advise is that -- if at all possible -- you should try to select a design process that by its nature generates a list of 2nd-order filters, rather than obfuscating the pole/zero positions by munging them all together into one Master Transfer Function which you then have to take apart into 2nd-order sections. Because, guaranteed, no matter what, sooner or later you will specify some bazzilionth-order filter, the accurate dissection of which will be doomed to failure due to numerical sensitivity of the coefficients. Of course, if you're using a filter design tool that won't let you break things apart, you just need to curse in a good natured sort of way, carry on, and hope that your phase response will be close enough. -- My liberal friends think I'm a conservative kook. My conservative friends think I'm a liberal kook. Why am I not happy that they have found common ground? Tim Wescott, Communications, Control, Circuits & Software http://www.wescottdesign.com
Reply by ●May 10, 20122012-05-10
>On Thu, 10 May 2012 12:40:12 -0400, robert bristow-johnson wrote: > >> On 5/10/12 12:23 PM, gretzteam wrote: >>> Hi, >>> I have an 8th order allpass filter and want to implement it as a >>> cascade of 2nd order ALLPASS sections. The usual 'split into SOS' >>> method doesn't keep the allpass property, >> >> turn the APF into LPF (ignoring the numerator), factor into SecondOrder>> Sections, then reflect each pole explicitly to the APF zero. >> >>> so I now have to manually look at the pole/zeros of the high order >>> filter, pick the pairs that are 2nd order friendly, and rebuild the >>> filter. >> >> do the same (rebuild the filter), but if you are guaranteed that it's >> APF to start with, knowledge of the zeros doesn't do anything for you. >> just factor the poles and reflect each pole to its corresponding zero, >> rather than trying to find the zero that matches the pole. >> >> >> >>> Is there a better way to do this - maybe this factoring with allpass >>> constraint has a 'name'. > >If there's anything that I'd add to Robert's advise, it is this: > >Your problem most likely stems from trying to factor high-degree >polynomials. Polynomial roots are notoriously sensitive to coefficient >variations, and get ever more so as the order gets higher. > >Yet IIR filters are sensitive to the _roots_ of their numerator and >denominator -- not the actual coefficients of the overall transfer >function. > >Which means that your problem with your gazzilionth-order numerator and >gazzilionth-order denominator isn't that you are getting your desired >pole positions but not your desired zeros -- it means that _both_ your >poles and your zeros are screwed up, but they just happen to be screwed >up differently. > >Using Robert's method is essentially saying "OK, I'm going to live with >the wrong pole positions (and hence a phase response that isn't to my >original design) as long as my amplitude stays flat. > >So _my_ advise is that -- if at all possible -- you should try to select >a design process that by its nature generates a list of 2nd-order >filters, rather than obfuscating the pole/zero positions by munging them >all together into one Master Transfer Function which you then have to >take apart into 2nd-order sections. Because, guaranteed, no matter what,>sooner or later you will specify some bazzilionth-order filter, the >accurate dissection of which will be doomed to failure due to numerical >sensitivity of the coefficients. > >Of course, if you're using a filter design tool that won't let you break >things apart, you just need to curse in a good natured sort of way, carry>on, and hope that your phase response will be close enough. > >-- >My liberal friends think I'm a conservative kook. >My conservative friends think I'm a liberal kook. >Why am I not happy that they have found common ground? > >Tim Wescott, Communications, Control, Circuits & Software >http://www.wescottdesign.com >Thanks for the comment. My filter design method cannot generate 2nd order section directly, but i'm only interested in 4th and 8th order filters, and so far Robert's method works nicely - and I went down from 23 lines of code + eyeballing to 4 simple lines!
Reply by ●May 10, 20122012-05-10
On Thu, 10 May 2012 14:22:23 -0500, gretzteam wrote:>>On Thu, 10 May 2012 12:40:12 -0400, robert bristow-johnson wrote: >> >>> On 5/10/12 12:23 PM, gretzteam wrote: >>>> Hi, >>>> I have an 8th order allpass filter and want to implement it as a >>>> cascade of 2nd order ALLPASS sections. The usual 'split into SOS' >>>> method doesn't keep the allpass property, >>> >>> turn the APF into LPF (ignoring the numerator), factor into Second > Order >>> Sections, then reflect each pole explicitly to the APF zero. >>> >>>> so I now have to manually look at the pole/zeros of the high order >>>> filter, pick the pairs that are 2nd order friendly, and rebuild the >>>> filter. >>> >>> do the same (rebuild the filter), but if you are guaranteed that it's >>> APF to start with, knowledge of the zeros doesn't do anything for you. >>> just factor the poles and reflect each pole to its corresponding zero, >>> rather than trying to find the zero that matches the pole. >>> >>> >>> >>>> Is there a better way to do this - maybe this factoring with allpass >>>> constraint has a 'name'. >> >>If there's anything that I'd add to Robert's advise, it is this: >> >>Your problem most likely stems from trying to factor high-degree >>polynomials. Polynomial roots are notoriously sensitive to coefficient >>variations, and get ever more so as the order gets higher. >> >>Yet IIR filters are sensitive to the _roots_ of their numerator and >>denominator -- not the actual coefficients of the overall transfer >>function. >> >>Which means that your problem with your gazzilionth-order numerator and >>gazzilionth-order denominator isn't that you are getting your desired >>pole positions but not your desired zeros -- it means that _both_ your >>poles and your zeros are screwed up, but they just happen to be screwed >>up differently. >> >>Using Robert's method is essentially saying "OK, I'm going to live with >>the wrong pole positions (and hence a phase response that isn't to my >>original design) as long as my amplitude stays flat. >> >>So _my_ advise is that -- if at all possible -- you should try to select >>a design process that by its nature generates a list of 2nd-order >>filters, rather than obfuscating the pole/zero positions by munging them >>all together into one Master Transfer Function which you then have to >>take apart into 2nd-order sections. Because, guaranteed, no matter >>what, > >>sooner or later you will specify some bazzilionth-order filter, the >>accurate dissection of which will be doomed to failure due to numerical >>sensitivity of the coefficients. >> >>Of course, if you're using a filter design tool that won't let you break >>things apart, you just need to curse in a good natured sort of way, >>carry > >>on, and hope that your phase response will be close enough. >> >>-- >>My liberal friends think I'm a conservative kook. My conservative >>friends think I'm a liberal kook. Why am I not happy that they have >>found common ground? >> >>Tim Wescott, Communications, Control, Circuits & Software >>http://www.wescottdesign.com >> >> > Thanks for the comment. My filter design method cannot generate 2nd > order section directly, but i'm only interested in 4th and 8th order > filters, and so far Robert's method works nicely - and I went down from > 23 lines of code + eyeballing to 4 simple lines!8th order is a stretch, but if you're getting sufficient precision to keep the pole damping ratios from changing too much, then you're probably golden. -- My liberal friends think I'm a conservative kook. My conservative friends think I'm a liberal kook. Why am I not happy that they have found common ground? Tim Wescott, Communications, Control, Circuits & Software http://www.wescottdesign.com
Reply by ●May 10, 20122012-05-10
On 5/10/12 4:11 PM, Tim Wescott wrote:> On Thu, 10 May 2012 14:22:23 -0500, gretzteam wrote: > >> Thanks for the comment. My filter design method cannot generate 2nd >> order section directly, but i'm only interested in 4th and 8th order >> filters, and so far Robert's method works nicely - and I went down from >> 23 lines of code + eyeballing to 4 simple lines! > > 8th order is a stretch,double-precision in MATLAB?? (it's gotta be something line MATLAB if gretz can do this in 4 lines of code.) factoring an 8th-order polynomial should not be a problem for floating-point with 53 bits of mantissa. or even stretching it. bump that up to a "bazillionth order" and any finite-precision arithmetic will fall short.> but if you're getting sufficient precision to > keep the pole damping ratios from changing too much, then you're probably > golden.sounds like gretz is doing something with Linkwitz-Riley. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●May 10, 20122012-05-10
On Thu, 10 May 2012 16:20:59 -0400, robert bristow-johnson wrote:> On 5/10/12 4:11 PM, Tim Wescott wrote: >> On Thu, 10 May 2012 14:22:23 -0500, gretzteam wrote: >> >>> Thanks for the comment. My filter design method cannot generate 2nd >>> order section directly, but i'm only interested in 4th and 8th order >>> filters, and so far Robert's method works nicely - and I went down >>> from 23 lines of code + eyeballing to 4 simple lines! >> >> 8th order is a stretch, > > double-precision in MATLAB?? (it's gotta be something line MATLAB if > gretz can do this in 4 lines of code.) factoring an 8th-order > polynomial should not be a problem for floating-point with 53 bits of > mantissa. or even stretching it. bump that up to a "bazillionth order" > and any finite-precision arithmetic will fall short.Well, probably not. But even an 8th-order will have problems if the Q is too high or if the normalized center frequency is low. An 8-th order polynomials with all the root magnitudes on the order of |z|=1 means that the factorization is going to (roughly) good to an error magnitude of 0.01. That would work out to a 10% frequency change in a pole at 700Hz, or a change in 0.1 in the damping ratio, if you were sampling at 44000.> >> but if you're getting sufficient precision to keep the pole damping >> ratios from changing too much, then you're probably golden. > > sounds like gretz is doing something with Linkwitz-Riley.-- My liberal friends think I'm a conservative kook. My conservative friends think I'm a liberal kook. Why am I not happy that they have found common ground? Tim Wescott, Communications, Control, Circuits & Software http://www.wescottdesign.com
Reply by ●May 10, 20122012-05-10
On 5/10/12 5:07 PM, Tim Wescott wrote:> On Thu, 10 May 2012 16:20:59 -0400, robert bristow-johnson wrote: > >> On 5/10/12 4:11 PM, Tim Wescott wrote: >>> On Thu, 10 May 2012 14:22:23 -0500, gretzteam wrote: >>> >>>> Thanks for the comment. My filter design method cannot generate 2nd >>>> order section directly, but i'm only interested in 4th and 8th order >>>> filters, and so far Robert's method works nicely - and I went down >>>> from 23 lines of code + eyeballing to 4 simple lines! >>> >>> 8th order is a stretch, >> >> double-precision in MATLAB?? (it's gotta be something line MATLAB if >> gretz can do this in 4 lines of code.) factoring an 8th-order >> polynomial should not be a problem for floating-point with 53 bits of >> mantissa. or even stretching it. bump that up to a "bazillionth order" >> and any finite-precision arithmetic will fall short. > > Well, probably not. But even an 8th-order will have problems if the Q is > too high or if the normalized center frequency is low. An 8-th order > polynomials with all the root magnitudes on the order of |z|=1 means that > the factorization is going to (roughly) good to an error magnitude of > 0.01. That would work out to a 10% frequency change in a pole at 700Hz, > or a change in 0.1 in the damping ratio, if you were sampling at 44000. >so you're saying (and you can well be right) that for an 8th-order factorization problem, that even with 53 meaningful bits of in words of any parameters and intermediate variables, that if the coefficients are nasty enough that there are zeros real close to |z|=1 (do they have to be close to each other for this to get pathological?), that you might expect only 7 or 8 meaningful bits in the real and imaginary parts of the roots? can't you precisely check the root by just plugging it into the original polynomial? oh, i get it: 0.01^8 is about 2^(-53). i guess this is why there is an Add-With-Carry instruction. i dunno if there is a MATLAB toolbox that does this, but i could imagine that there is some programs out there that can factor nasty polynomials by making really long fixed or floating-point words. but i don't know how you would write the code in C, because you need Add-With-Carry. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●May 11, 20122012-05-11
>On 5/10/12 5:07 PM, Tim Wescott wrote: >> On Thu, 10 May 2012 16:20:59 -0400, robert bristow-johnson wrote: >> >>> On 5/10/12 4:11 PM, Tim Wescott wrote: >>>> On Thu, 10 May 2012 14:22:23 -0500, gretzteam wrote: >>>> >>>>> Thanks for the comment. My filter design method cannot generate 2nd >>>>> order section directly, but i'm only interested in 4th and 8th order >>>>> filters, and so far Robert's method works nicely - and I went down >>>>> from 23 lines of code + eyeballing to 4 simple lines! >>>> >>>> 8th order is a stretch, >>> >>> double-precision in MATLAB?? (it's gotta be something line MATLAB if >>> gretz can do this in 4 lines of code.) factoring an 8th-order >>> polynomial should not be a problem for floating-point with 53 bits of >>> mantissa. or even stretching it. bump that up to a "bazillionthorder">>> and any finite-precision arithmetic will fall short. >>Those are phase compensation filters. I didn't see any problems in quantizing the results and keeping the phase acceptable, although I've only tried a 4th order so far. Actually, I'm quite surprised by the robustness of the quantization - of course I'm using a miminum-multiplier allpass structure and not a direct-form. This way it stays allpass no matter the coefficient quantization - kind of the dual of what rbj proposed for the factoring into SOS, but on the implementation side.
