How can I scale the coefficients of an second order section IIR filter so they are within -1 to +0.99999988079071 ? When I calculate a filter in MATLAB fdatool I get some coefficients that are > 1.
How to scale second order IIR filter koefficients?
Started by ●August 15, 2006
Reply by ●August 15, 20062006-08-15
"Gery" <gery@ddd.com> writes:> How can I scale the coefficients of an second order section IIR filter so > they are within -1 to +0.99999988079071 ? > > > > When I calculate a filter in MATLAB fdatool I get some coefficients that are > > 1.Ideally, you can't, because the "1" in the denominator of a SOS, expressed as H(z)= (b2*z^2 + b1*z + b0) / (1 + a1*z + a0), must always be 1. If you're asking how to represent the b's and a's in fixed-point form, then that's a different questtion. -- % Randy Yates % "Bird, on the wing, %% Fuquay-Varina, NC % goes floating by %%% 919-577-9882 % but there's a teardrop in his eye..." %%%% <yates@ieee.org> % 'One Summer Dream', *Face The Music*, ELO http://home.earthlink.net/~yatescr
Reply by ●August 15, 20062006-08-15
Randy Yates wrote:> "Gery" <gery@ddd.com> writes: > > >>How can I scale the coefficients of an second order section IIR filter so >>they are within -1 to +0.99999988079071 ? >> >> >> >>When I calculate a filter in MATLAB fdatool I get some coefficients that are >> > 1. > > > Ideally, you can't, because the "1" in the denominator of a SOS, > expressed as > > H(z)= (b2*z^2 + b1*z + b0) / (1 + a1*z + a0), > > must always be 1. > > If you're asking how to represent the b's and a's in fixed-point > form, then that's a different questtion.I think you mean (1*z^2 + a1*z + a0). -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Posting from Google? See http://cfaj.freeshell.org/google/ "Applied Control Theory for Embedded Systems" came out in April. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●August 15, 20062006-08-15
Gery wrote:> How can I scale the coefficients of an second order section IIR filter so > they are within -1 to +0.99999988079071 ? > > > > When I calculate a filter in MATLAB fdatool I get some coefficients that are > > 1. > >Folks still use the 56xxx? Cool! So you have a transfer function that goes like: b2*z^2 + b1*z + b0 H(z) = ------------------, z^2 + a1*z + a0 and a1 is greater than 1, right? If you turn this into the difference equation x[0] = -a1*x[1] - a0*x[2] + b2*y[0] + b1*y[1] + b0*y[2] Presumably you're doing this in fractional arithmetic, which is why you have the concern about coefficients, and it's a1 that's giving you the trouble. You can do several things to fix this. Probably the easiest thing to do is to realize your difference equation as: x[0] = x[1] - (a1 + 1)*x[1] - a0*x[2] + b2*y[0] + b1*y[1] + b0*y[2]. Since your a1 is generally going to be in the range (-2, -1) this works nicely. The second easiest thing to do is realize it as: x[0] = 2 * (-a1/2*x[1]) - a0*x[2] + b2*y[0] + b1*y[1] + b0*y[2]. The multiply by 2 can be done as a shift at the cost of one machine cycle, and you're not losing precision on a0, which is where it really matters. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Posting from Google? See http://cfaj.freeshell.org/google/ "Applied Control Theory for Embedded Systems" came out in April. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●August 15, 20062006-08-15
Randy Yates wrote:>>How can I scale the coefficients of an second order section IIR filter so >>they are within -1 to +0.99999988079071 ? >> > > Ideally, you can't, because the "1" in the denominator of a SOS, > expressed as > > H(z)= (b2*z^2 + b1*z + b0) / (1 + a1*z + a0), > > must always be 1. > > If you're asking how to represent the b's and a's in fixed-point > form, then that's a different questtion.In the 2nd order IIR denominator, b1 is usually at the order of 2. The easiest thing to do if you have to use the fractional math is to repeat the MAC operation with b1 twice. This ensures b1 is less then +/- 1. Vladimir Vassilevsky DSP and Mixed Signal Design Consultant http://www.abvolt.com
Reply by ●August 15, 20062006-08-15
Tim Wescott wrote:> Gery wrote: > > > How can I scale the coefficients of an second order section IIR filter so > > they are within -1 to +0.99999988079071 ? > > > > > > > > When I calculate a filter in MATLAB fdatool I get some coefficients that are > > > 1. > > > > > Folks still use the 56xxx? Cool!digidesign Pro Tools TDM still uses DSP563xx. audio folks have been speculating for years how long it would be until they switch to SHArC or something else. but i didn't infer from the OP it was 56K.> > So you have a transfer function that goes like: > > b2*z^2 + b1*z + b0 > H(z) = ------------------, > z^2 + a1*z + a0 > > and a1 is greater than 1, right? > > If you turn this into the difference equation > > x[0] = -a1*x[1] - a0*x[2] + b2*y[0] + b1*y[1] + b0*y[2] > > Presumably you're doing this in fractional arithmetic, which is why you > have the concern about coefficients, and it's a1 that's giving you the > trouble. You can do several things to fix this. Probably the easiest > thing to do is to realize your difference equation as: > > x[0] = x[1] - (a1 + 1)*x[1] - a0*x[2] + b2*y[0] + b1*y[1] + b0*y[2]. > > Since your a1 is generally going to be in the range (-2, -1) this works > nicely.generally -2 < a1 = 2*|p|*cos(theta) < +2 and 0 <= a2 = |p|^2 < 1 so usually when i see coefficients that look like "+0.99999988079071" i usually think that it's "a2" for a very high Q filters.> The second easiest thing to do is realize it as: > > x[0] = 2 * (-a1/2*x[1]) - a0*x[2] + b2*y[0] + b1*y[1] + b0*y[2].i think you got your inputs and outputs symbols switched, Tim. i'm switching some other symbols: b0 + b1*z^-1 + b2*z^-2 Y(z)/X(z) = H(z) = ------------------------- 1 + a1*z^-1 + a2^-2 because the feedforward coefs can be practically anything, i usually recommend: y[n] = 2*( -a1/2*y[n-1]) - a2/2*y[n-2] + N/2*(b0/N*x[n] + b1/N*x[n-1] + a2/N*y[n-2]) ) where N is another power of 2, usually 4 or 8 or maybe 16. those the coefs stored is as shown and the multiplying by 2 or N/2 is done by a couple of ASL instructions operating on the double wide accumulator. r b-j
Reply by ●August 15, 20062006-08-15
Tim Wescott <tim@seemywebsite.com> writes:> Randy Yates wrote: > >> "Gery" <gery@ddd.com> writes: >> >>> How can I scale the coefficients of an second order section IIR >>> filter so they are within -1 to +0.99999988079071 ? >>> >>> >>> >>> When I calculate a filter in MATLAB fdatool I get some coefficients >>> that are > 1. >> Ideally, you can't, because the "1" in the denominator of a SOS, >> expressed as >> H(z)= (b2*z^2 + b1*z + b0) / (1 + a1*z + a0), >> must always be 1. If you're asking how to represent the b's and a's >> in fixed-point >> form, then that's a different questtion. > > I think you mean (1*z^2 + a1*z + a0).Yes, that is what I meant. Thank you for the correction, Tim. -- % Randy Yates % "The dreamer, the unwoken fool - %% Fuquay-Varina, NC % in dreams, no pain will kiss the brow..." %%% 919-577-9882 % %%%% <yates@ieee.org> % 'Eldorado Overture', *Eldorado*, ELO http://home.earthlink.net/~yatescr
Reply by ●August 15, 20062006-08-15
robert bristow-johnson wrote:> Tim Wescott wrote: > >>Gery wrote: >> >> >>>How can I scale the coefficients of an second order section IIR filter so >>>they are within -1 to +0.99999988079071 ? >>> >>> >>> >>>When I calculate a filter in MATLAB fdatool I get some coefficients that are >>> > 1. >>> >>> >> >>Folks still use the 56xxx? Cool! > > > digidesign Pro Tools TDM still uses DSP563xx. audio folks have been > speculating for years how long it would be until they switch to SHArC > or something else. > > but i didn't infer from the OP it was 56K. >0.99999988079071 = 1 - 2^-23, from which I inferred a 24-bit word length, from which I inferred a 56K. It's a weak chain of reasoning, I'll admit -- but are there any other 24-bit DSP chips out there? -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Posting from Google? See http://cfaj.freeshell.org/google/ "Applied Control Theory for Embedded Systems" came out in April. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●August 15, 20062006-08-15
Vladimir Vassilevsky wrote: ...> In the 2nd order IIR denominator, b1 is usually at the order of 2. The > easiest thing to do if you have to use the fractional math is to repeat > the MAC operation with b1 twice. This ensures b1 is less then +/- 1.Vladimir's native language probably isn't English. I don't fault him, but rather the native anglophones his sentence was modeled on. The MAC operation is _done_ twice; that means it is repeated _once_. Vladimir has gotten a high (but silent) mark from me in the past for writing "there is an alternative" rather than the all-too-common "there is another alternative" when "there is another way" is meant. I salute him now. May he join the ranks of Joseph Conrad and Vladimir Nabokov in the ESL Hall Of Fame. :-) Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●August 15, 20062006-08-15
Tim Wescott wrote:> robert bristow-johnson wrote: > > > Tim Wescott wrote: > > > >>Gery wrote: > >> > >> > >>>How can I scale the coefficients of an second order section IIR filter so > >>>they are within -1 to +0.99999988079071 ? > >>> > >>> > >>> > >>>When I calculate a filter in MATLAB fdatool I get some coefficients that are > >>> > 1. > >>> > >> > >>Folks still use the 56xxx? Cool! > > > > > > digidesign Pro Tools TDM still uses DSP563xx. audio folks have been > > speculating for years how long it would be until they switch to SHArC > > or something else. > > > > but i didn't infer from the OP it was 56K. > > > 0.99999988079071 = 1 - 2^-23, from which I inferred a 24-bit word > length, from which I inferred a 56K. It's a weak chain of reasoning, > I'll admit -- but are there any other 24-bit DSP chips out there?not common. there are some cores, and (i didn't know this before) there are even 56K clone cores: http://www.chipestimate.com/ip.php?id=5812 (ya gotta register, big headache). but there are others like http://www.clarkspur.com/sum2480a.htm . one thing that i still lament is that Motorola had the audio DSP chip market locked and they dropped the ball. they needed to do something to response to ADI especially when the ADSP-21065L came out. now, so much audio processing is done native that it really doesn't matter. r b-j
