| Author |
Message |
Fred Nach
Guest
|
 Posted:
Wed Nov 30, 2005 5:16 pm Post subject:
Fixed Point IIR implementation |
|
|
Hi pals,
I would like to implement a IIR Biquad filter using the fixed point
arithmetics...
Hence to reduced the intermediate states I plan to use the following trick:
s(k) = x(k) -a1*s(k-1) -a2*s(k-2)
y(n) = b0*s(k) + b1*s(k-1) + b2*s(k-2)
The I can compute each y, saving only 2 states (s(k-1) aénd s(k-2))...
BUT ... in order to scale my coefficients or input, I need to know what
are the boundaries of the s(k) serie... in order to find a proper Fixed
Point format for my coefficients..
So my questions are :
* for which conditions, s(k) is bounded
* If the s(k) is bounded, what are the boundaries?
Thanks in advance for your help (and sorry for my poor englsih)
Fred. |
|
| Back to top |
|
 |
Jerry Avins
Guest
|
| Posted:
Wed Nov 30, 2005 5:17 pm Post subject:
Re: Fixed Point IIR implementation |
|
|
Fred Nach wrote:
| Quote: | Hi pals,
I would like to implement a IIR Biquad filter using the fixed point
arithmetics...
Hence to reduced the intermediate states I plan to use the following trick:
s(k) = x(k) -a1*s(k-1) -a2*s(k-2)
y(n) = b0*s(k) + b1*s(k-1) + b2*s(k-2)
The I can compute each y, saving only 2 states (s(k-1) aénd s(k-2))...
BUT ... in order to scale my coefficients or input, I need to know what
are the boundaries of the s(k) serie... in order to find a proper Fixed
Point format for my coefficients..
So my questions are :
* for which conditions, s(k) is bounded
* If the s(k) is bounded, what are the boundaries?
Thanks in advance for your help (and sorry for my poor englsih)
|
Your English is good enough, I think. The s[k] are data you supply. They
can't be larger or smaller than the representation allows, but they may
have smaller limits.
I see no feedback in the formulas you give. Is that correct?
Jerry
--
Engineering is the art of making what you want from things you can get.
ŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻ |
|
| Back to top |
|
|
robert bristow-johnson
Guest
|
| Posted:
Wed Nov 30, 2005 11:57 pm Post subject:
Re: Fixed Point IIR implementation |
|
|
Fred Nach wrote:
| Quote: |
I would like to implement a IIR Biquad filter using the fixed point
arithmetics...
Hence to reduced the intermediate states I plan to use the following trick:
s(k) = x(k) -a1*s(k-1) -a2*s(k-2)
y(n) = b0*s(k) + b1*s(k-1) + b2*s(k-2)
The I can compute each y, saving only 2 states (s(k-1) aénd s(k-2))...
|
this is the Direct 2 Form (sometimes called the "Direct II Canonical
Form"). unrecommeded for floating-point (if the Q is quite high, you
end up subtracting numbers very close to each other and losing
precision) and *highly* unrecommended for fixed-point (saturation
clipping will be your lot).
try the Direct 1 Form:
y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] - a2*y[n-2]
every product on the right side of the = sign is a double precision
fixed point and all of those products should be added together in
double precision. this is trivial in the Mot 56K and the ADI SHArC and
maybe in the TI fixed-pointers.
| Quote: | BUT ... in order to scale my coefficients or input, I need to know what
are the boundaries of the s(k) serie... in order to find a proper Fixed
Point format for my coefficients..
|
your coefs are defined strictly in terms of the frequency response you
want. the gross range of a1 is -2 to +2 and a2 is -1 to +1 for any
stable biquad filter. the b0, b1, b2 coefs can be nearly anything but
you will have to choose a range and possibly do some scaling using the
arithmetic shift operation on the result before saving to y[n].
| Quote: | So my questions are :
* for which conditions, s(k) is bounded
* If the s(k) is bounded, what are the boundaries?
|
that is a purely artificial construct. are your fixed-point signal
values considered 16 bit signed integers? then it's -32768 <= x[n] <
+32768 . are they normalized? then it's -1 <= x[n]< +1. their range
is whatever you want them to be.
| Quote: | Thanks in advance for your help (and sorry for my poor englsih)
|
as Jerry said, your English is fine.
r b-j |
|
| Back to top |
|
|
Noway2
Guest
|
| Posted:
Thu Dec 01, 2005 12:01 am Post subject:
Re: Fixed Point IIR implementation |
|
|
| Quote: | I see no feedback in the formulas you give. Is that correct?
|
The aX coeffecients should be the feedback terms and the bX
coefficients should be the feed forward terms (from drawing a direct
form II simulation diagram). The s(k) terms are the inputs to the
delay buffers, which will be a function of the input, feed forward and
feedback terms.
I see the question he is asking, how to determine if these terms
overflow the fixed point number format in use, and how to determine a
proper scaling value to keep it from happening. Unfortunately, I don't
know how to answer this question, but in the project I am presently
working on I have encountered the same problem, so I am hoping someone
has a good answer. I got around the problem (? maybe ?) by using a
direct form I structure.
I have seen Matlab scripts that run an impulse response till it
stabilizes and then sum up the impulse responses and use the maximum
value as a scale factor, though I don't know if this is a correct
solution. |
|
| Back to top |
|
|
Jerry Avins
Guest
|
| Posted:
Thu Dec 01, 2005 12:36 am Post subject:
Re: Fixed Point IIR implementation |
|
|
Noway2 wrote:
| Quote: | I see no feedback in the formulas you give. Is that correct?
The aX coeffecients should be the feedback terms and the bX
coefficients should be the feed forward terms (from drawing a direct
form II simulation diagram). The s(k) terms are the inputs to the
delay buffers, which will be a function of the input, feed forward and
feedback terms.
I see the question he is asking, how to determine if these terms
overflow the fixed point number format in use, and how to determine a
proper scaling value to keep it from happening. Unfortunately, I don't
know how to answer this question, but in the project I am presently
working on I have encountered the same problem, so I am hoping someone
has a good answer. I got around the problem (? maybe ?) by using a
direct form I structure.
I have seen Matlab scripts that run an impulse response till it
stabilizes and then sum up the impulse responses and use the maximum
value as a scale factor, though I don't know if this is a correct
solution.
|
I expect that the b's are the denominator terms. Input samples are s[n],
output and feedback samples are y[n]. I see no coefficients operating on
any y[n].
I see what he's asking too. His problem lies with internal states, and I
don't know how to deal with that in the abstract. If there's a general
case, I haven't seen it.
Jerry
--
Engineering is the art of making what you want from things you can get.
ŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻ |
|
| Back to top |
|
|
Andrew Reilly
Guest
|
| Posted:
Thu Dec 01, 2005 1:17 am Post subject:
Re: Fixed Point IIR implementation |
|
|
On Wed, 30 Nov 2005 11:01:31 -0800, Noway2 wrote:
| Quote: | I have seen Matlab scripts that run an impulse response till it
stabilizes and then sum up the impulse responses and use the maximum
value as a scale factor, though I don't know if this is a correct
solution.
|
It's pretty close. This is the main reason why these "direct form 2"
implmentations are not preferred. A direct form 1 implmentation will use
more memory but (a) memory is relatively cheap and (b) the values stored
will be as constrained as the input and output values.
--
Andrew |
|
| Back to top |
|
|
robert bristow-johnson
Guest
|
| Posted:
Thu Dec 01, 2005 1:17 am Post subject:
Re: Fixed Point IIR implementation |
|
|
Jerry Avins wrote:
| Quote: |
I expect that the b's are the denominator terms.
|
no, Jerry, the a's are in the denominator. they reversed that
convention a while back. in the 70's it was a_k on top and b_k on
bottom, but, i think because they think that poles are more important,
most texts have it switched around now.
r b-j |
|
| Back to top |
|
|
robert bristow-johnson
Guest
|
| Posted:
Thu Dec 01, 2005 1:17 am Post subject:
Re: Fixed Point IIR implementation |
|
|
Jerry Avins wrote:
| Quote: | robert bristow-johnson wrote:
Jerry Avins wrote:
I expect that the b's are the denominator terms.
no, Jerry, the a's are in the denominator.
I knew that too, but old habits sometimes bubble up. The important point
was that Fred N. has no y terms to the right of an equal sign, so I
don't see any recursion.
|
there is recursion with the intermediate s[k] signal:
Fred Nach wrote:
| Quote: | I would like to implement a IIR Biquad filter using the fixed point
arithmetics...
Hence to reduced the intermediate states I plan to use the following trick:
s(k) = x(k) -a1*s(k-1) -a2*s(k-2)
y(n) = b0*s(k) + b1*s(k-1) + b2*s(k-2)
|
this is the canonical Direct Form 2. poles before zeros. problem is,
if there is a high Q filter, those poles will amplify the signal so
much that it will saturate s[k] before the zeros get to beat the signal
back down. DF2 ain't particularly good, especially for fixed-point
arithmetic.
r b-j |
|
| Back to top |
|
|
Jerry Avins
Guest
|
| Posted:
Thu Dec 01, 2005 1:17 am Post subject:
Re: Fixed Point IIR implementation |
|
|
robert bristow-johnson wrote:
| Quote: | Jerry Avins wrote:
I expect that the b's are the denominator terms.
no, Jerry, the a's are in the denominator. they reversed that
convention a while back. in the 70's it was a_k on top and b_k on
bottom, but, i think because they think that poles are more important,
most texts have it switched around now.
r b-j
|
I knew that too, but old habits sometimes bubble up. The important point
was that Fred N. has no y terms to the right of an equal sign, so I
don't see any recursion. Somewhere along the line I wrote that. Twice, I
think.
Jerry
--
Engineering is the art of making what you want from things you can get.
ŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻ |
|
| Back to top |
|
|
Jerry Avins
Guest
|
| Posted:
Thu Dec 01, 2005 9:16 am Post subject:
Re: Fixed Point IIR implementation |
|
|
robert bristow-johnson wrote:
| Quote: | Jerry Avins wrote:
robert bristow-johnson wrote:
Jerry Avins wrote:
I expect that the b's are the denominator terms.
no, Jerry, the a's are in the denominator.
I knew that too, but old habits sometimes bubble up. The important point
was that Fred N. has no y terms to the right of an equal sign, so I
don't see any recursion.
there is recursion with the intermediate s[k] signal:
|
OK. Back to the books with my thinking cap on.
Jerry
--
Engineering is the art of making what you want from things you can get.
ŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻ |
|
| Back to top |
|
|
robert bristow-johnson
Guest
|
| Posted:
Thu Dec 01, 2005 9:16 am Post subject:
Re: Fixed Point IIR implementation |
|
|
Jerry Avins wrote:
| Quote: | robert bristow-johnson wrote:
Jerry Avins wrote:
....
I knew that too, but old habits sometimes bubble up. The important point
was that Fred N. has no y terms to the right of an equal sign, so I
don't see any recursion.
there is recursion with the intermediate s[k] signal:
OK. Back to the books with my thinking cap on.
|
no need. here's a simple way to look at it:
s[k] = x[k] - a1*s[k-1] - a2*s[k-2]
y[k] = b0*s[k] + b1*s[k-1] + b2*s[k-2]
think of this as two filters in cascade. the first with input x[k] and
output s[k]. the second with input s[k] and output y[k].
the first is a all-pole filter (actually there are 2 zeros at the
origin, but don't worry about them) and the second is a all zero
filter, an FIR (like an FIR, there are 2 poles at the origin, but the
other 2 zeros kill 'em). the transfer function of the first:
S(z)/X(z) = 1/(1 + a1*z^-1 + a2*z-2)
the transfer function of the second
Y(z)/S(z) = b0 + b1*z^-1 + b2*z^-2
and the transfer function of the whole sheee-bang is
H(z) = Y(z)/X(z) = ( Y(z)/S(z) ) * ( S(z)/X(z) )
and the reason it has only two states instead of four is that the
states for the first are identical to the states of the second (both
delays on s[k]) so they can be combined into a single delay line.
r b-j |
|
| Back to top |
|
|
Fred Nach
Guest
|
| Posted:
Thu Dec 01, 2005 5:16 pm Post subject:
Re: Fixed Point IIR implementation |
|
|
Me again ...
Ok, so i've just launch a Matlab plot of the filter derived from the
s(k) sequence :
s(k) = x(k) -a1*s(k-1) -a2*s(k-2)
hence considering x(k) as input and s(k) as output
H(Z) = 1/(1 + a1*z^-1 + a2*z^-2)
the frequency response give a curve that decrease gradualy, *BUT* start
witha gain of about 80 dB ... hence it means (as long as i
understand...) that for almost constant x(k) (=low freq) the gain is
very important so s(k) will take big values ... (gain of 80dB is like
multiplying by 10 000 !) hence it for a Fixed Point architecture it
requires ceil(log2(10000)) = 14 guard bits !!!
So yes... it's definitively not a good solution to implement a filter
like that ... and I will do it using the direct form 1.
Salut !
Fred |
|
| Back to top |
|
|
robert bristow-johnson
Guest
|
| Posted:
Thu Dec 01, 2005 5:16 pm Post subject:
Re: Fixed Point IIR implementation |
|
|
Fred Nach wrote:
| Quote: |
It seems that my previous post was not clear enough
|
it was clear enough for me.
....
| Quote: | * The filter is indeed stable hence |a1| <= 2 and |a2| <= 1
I think I understood that this implementation is not recommended because
it compute first the feedback terms and leads to big variations...
|
it means that s[n] is too large to fit in you 16 bit words. if you
represent them as 32 bit words, then you have to do double precision
multiplies when you multiply be b0, b1, b2.
| Quote: | But is there a way to retrieve the boudary of s[n] mathematicaly ?
I think that as long as a1 and a2 are given for a stable filter s[n]
should be bounded...
|
it is. the worst case is that if a1 and a2 are known, you can obtain
the impulse response for the all-pole filter
H_s(z) = 1/(1 + a1*z^-1 + a2*z^-2)
call that impulse response h_s[n].
the worst case for the maximum value for |s[n]| is what happens if x[n]
shares the same sign as s[n] (for each value of n) and is at it at its
maximum magnitude. from the convolution summation
s[n] = SUM{ h_s[k] * x[n-k] }
n
max |s[n]| = max |x[n]| * SUM{ |h_s[n]| }
n n n (max's and SUM over all n)
you might find that max gain, SUM{ |h_s[n]| }, to be a very large
value. much larger than 1. if you were to reduce you input to such a
degree that you guarantee that s[n] never saturates, you will likely
reduce
| Quote: | For my program, this implementation is very interesting cause i code in
Assembly and It save a lot of MIPS to only read 2 states variable (s1
and s2) then 4 (x1 x2 y1 y2)...
|
there is no reason for the DF1 to require more MIPS in the inner loop
than the DF2. at least for cascaded biquads. if you're doing only one
biquad, that is your overall filter order is 2, then it will cost a
little more (but the cost is worth it).
the DF1 requires 2 additional states but, if you cascade 2nd order
biquads, the two output states can be shared with the 2 input states of
the next section. so if you have N sections(an order 2N filter), the
number of states you need for DF1 is 2N+2 only 2 more than DF2. not a
very high price to pay.
the DF2 is useless in your case (assuming that your filter will have
some decent Q making the poles very close to the unit circle). but, if
you like learning this yourself (we call this "learning the hard way"),
feel free to implement it. but you will find that either you
pre-scaling of the input will drive your signal into the noise floor
or,if you do not pre-scale the input, your intermediate state signal,
s[n], will saturate. only with double precision (that will cost more
MIPS than DF1) will you avoid that.
r b-j |
|
| Back to top |
|
|
robert bristow-johnson
Guest
|
| Posted:
Thu Dec 01, 2005 5:16 pm Post subject:
Re: Fixed Point IIR implementation |
|
|
Fred Nach wrote:
| Quote: | Me again ...
Ok, so i've just launch a Matlab plot of the filter derived from the
s(k) sequence :
s(k) = x(k) -a1*s(k-1) -a2*s(k-2)
hence considering x(k) as input and s(k) as output
H(Z) = 1/(1 + a1*z^-1 + a2*z^-2)
the frequency response give a curve that decrease gradualy, *BUT* start
witha gain of about 80 dB ... hence it means (as long as i
understand...) that for almost constant x(k) (=low freq) the gain is
very important so s(k) will take big values ... (gain of 80dB is like
multiplying by 10 000 !) hence it for a Fixed Point architecture it
requires ceil(log2(10000)) = 14 guard bits !!!
So yes... it's definitively not a good solution to implement a filter
like that ... and I will do it using the direct form 1.
|
our posts "crossed in the mail". i'm glad that you "get it" now. :-)
when using the DF1, make sure that all five terms (b0*x[n] b1*x[n-1],
etc) are added together in a double precision accumulator before
quantizing that back to a single precision output. also, consider
"noise shaping" or "fraction saving" (Google Groups the latter term or
look up the DC Blocking filter in fixed-point at dspguru.com). with
only 16 bits, you'll probably need it.
r b-j |
|
| Back to top |
|
|
Just Cocky
Guest
|
| Posted:
Thu Dec 01, 2005 5:16 pm Post subject:
Re: Fixed Point IIR implementation |
|
|
On 1 Dec 2005 08:35:02 -0800, "robert bristow-johnson"
<rbj@audioimagination.com> wrote:
| Quote: |
Fred Nach wrote:
Me again ...
Ok, so i've just launch a Matlab plot of the filter derived from the
s(k) sequence :
s(k) = x(k) -a1*s(k-1) -a2*s(k-2)
hence considering x(k) as input and s(k) as output
H(Z) = 1/(1 + a1*z^-1 + a2*z^-2)
the frequency response give a curve that decrease gradualy, *BUT* start
witha gain of about 80 dB ... hence it means (as long as i
understand...) that for almost constant x(k) (=low freq) the gain is
very important so s(k) will take big values ... (gain of 80dB is like
multiplying by 10 000 !) hence it for a Fixed Point architecture it
requires ceil(log2(10000)) = 14 guard bits !!!
So yes... it's definitively not a good solution to implement a filter
like that ... and I will do it using the direct form 1.
our posts "crossed in the mail". i'm glad that you "get it" now. :-)
when using the DF1, make sure that all five terms (b0*x[n] b1*x[n-1],
etc) are added together in a double precision accumulator before
quantizing that back to a single precision output. also, consider
"noise shaping" or "fraction saving" (Google Groups the latter term or
look up the DC Blocking filter in fixed-point at dspguru.com). with
only 16 bits, you'll probably need it.
|
Isn't this a case where a lattice filter implementaion helps? |
|
| Back to top |
|
|
|
|
|
|
FAQ
Memberlist
Register
Profile
Log in to check your private messages
Log in