Searching \ for '[PIC]: Software low pass filtering' in subject line. () Help us get a faster server
FAQ page: www.piclist.com/techref/microchip/math/filter.htm?key=filter
Search entire site for: 'Software low pass filtering'.

Exact match. Not showing close matches.
'[PIC]: Software low pass filtering'
2002\06\13@195542 by   Hello All,

some months ago Olin posted the formula that changed DSP world: :)
Filt <-- Filt + (New - Filt)*FFrac

Saddly, my 8 bit implementation have a problem that I can't surpass:
filt = filt + ((sig16)new - (sig16)filt) / 2;

The value never goes to 255 because of integer roundings:
filt = 254 + (255 - 254) / 2;
filt = 254 + (1)/2;
filt = 254;

There's a way to change this? I've tried to raise the resolution making:
filt = filt + ((sig16)new*2 - (sig16)filt*2 + 1) / 4;

But got same results because:
filt = 254 + ((sig16)255*2 - (sig16)254*2 + 1) / 4;
filt = 254 + (510 - 508 + 1) / 4;
filt = 254 + (510 - 508 + 1) / 4;
filt = 254 + (1) / 4;

Best regards,

Brusque

-----------------------------------------------------------------
Edson Brusque                 C.I.Tronics Lighting Designers Ltda
Research and Development               Blumenau  -  SC  -  Brazil
Say NO to HTML mail                          http://www.citronics.com.br
-----------------------------------------------------------------

--
http://www.piclist.com hint: The PICList is archived three different
ways.  See http://www.piclist.com/#archives for details.   Brusque --

> some months ago Olin posted the formula that changed DSP world: :)
>     Filt <-- Filt + (New - Filt)*FFrac
>
> Saddly, my 8 bit implementation have a problem that I can't surpass:
>     filt = filt + ((sig16)new - (sig16)filt) / 2;
>
> The value never goes to 255 because of integer roundings:
>     filt = 254 + (255 - 254) / 2;
>     filt = 254 + (1)/2;
>     filt = 254;
>
> There's a way to change this? I've tried to raise the resolution making:
>     filt = filt + ((sig16)new*2 - (sig16)filt*2 + 1) / 4;
>
> But got same results because:
>     filt = 254 + ((sig16)255*2 - (sig16)254*2 + 1) / 4;
>     filt = 254 + (510 - 508 + 1) / 4;
>     filt = 254 + (510 - 508 + 1) / 4;
>     filt = 254 + (1) / 4;

I think you correctly recognized that you need to round rather than
truncate the result, but you only rounded one term instead of the total:

Original:
filt = filt + (new - filt) / 2;

Rounded:
filt = ((filt*2 + (new*2 - filt*2) / 2) + 1) / 2;
filt = ((254*2 + (255*2 - 254*2) / 2) + 1) / 2;
filt = ((508 + (510 - 508) / 2) + 1) / 2;
filt = ((508 + 2 / 2) + 1) / 2;
filt = (509 + 1) / 2;
filt = 510 / 2;
filt = 255;

I hope this helps.

-- Dave Tweed

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body  Why do you care if it gets to 255 or just 254 ?

The error is pretty small.

Tal

{Quote hidden}

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body   On Thu, 13 Jun 2002, Edson Brusque (listas) wrote:

> Hello All,
>
>
>     some months ago Olin posted the formula that changed DSP world: :)
>         Filt <-- Filt + (New - Filt)*FFrac
>
>     Saddly, my 8 bit implementation have a problem that I can't surpass:
>         filt = filt + ((sig16)new - (sig16)filt) / 2;

The filter can be written like so:

Filt' = Filt + (New - Filt)*FFrac
= Filt(1-FFrac) + New*FFrac

Often, the FFrac is the ratio of two integers such that:

FFrac = B/(A+B)
1-FFrac = A/(A+B)

The sum of FFrac and (1-FFrac) is  (A+B)/(A+B) or 1.

Using this nomenclature, the filter can be rewritten:

Filt' = (A*Filt  + B*New)/(A+B)

So far nothing new, just a rearrangement of terms.

Now to apply rounding, there are few tricks that you can perform.

Filt' = (A*Filt + B*New + A*New)/(A+B) - A*New/(A+B)

Filt' = New + (Filt-New) * A/(A+B)        ........... (1)

>
>     The value never goes to 255 because of integer roundings:
>         filt = 254 + (255 - 254) / 2;
>         filt = 254 + (1)/2;
>         filt = 254;

In this filter, A=B=1. If you use equation 1:

filt = 255 + (254 - 255) / 2
= 255 - 1/2
= 255

Whoopee!

This rounding also works for 0: New=0, Filt=1

Filt' = 0 + (1-0)/2
= 1/2
= 0

One way to implement this would be:

; W = New

subwf  Filt,W    ;W = Filt-New
subwf  Filt,F    ;Swap New and Filtered values

movwf  temp      ;Save the difference
rlf    temp,w    ;put sign bit into carry
rrf    temp,w    ;signed right shift
btfsc  temp,7    ;If the difference is negative
addlw 1         ;  then round down
subwf  Filt,F

This implementation does have a subtle bug. If the difference between the
new and filtered values is greater than 0x7f, then the filter has a
spurious response. To fix this, you need to check the signs of the
difference versus the sign of the filt or new value. It can be done in 4
more instructions.

BTW, an efficient implementation of this class of recursive, one-pole
digital filters can be found:

http://www.dattalo.com/technical/software/pic/twist.asm

Scott

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body  On Thu, 13 Jun 2002, Edson Brusque (listas) wrote:

> Hello All,
>
>
>     some months ago Olin posted the formula that changed DSP world: :)
>         Filt <-- Filt + (New - Filt)*FFrac
>
>     Saddly, my 8 bit implementation have a problem that I can't surpass:
>         filt = filt + ((sig16)new - (sig16)filt) / 2;

This is because of the parity of 'filt' goes lost. I'd say filt should
inherit the parity of 'new'.

Regards,
Imre

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body  Hi,

Edson Brusque wrote:
>    some months ago Olin posted the formula that changed DSP world: :)
>        Filt <-- Filt + (New - Filt)*FFrac
>
>    Saddly, my 8 bit implementation have a problem that I can't
surpass:
>        filt = filt + ((sig16)new - (sig16)filt) / 2;
>
>    The value never goes to 255 because of integer roundings:
>        filt = 254 + (255 - 254) / 2;
>        filt = 254 + (1)/2;
>        filt = 254;
<snip>

Well this is one of the drawbacks with (standard) IIR filters, you're
only
affecting 1/x part of the difference to the output each pass.
In this case it easy:

Change:        Filt <-- Filt + (New - Filt)*FFrac

to:         Filt <-- Filt + (New - Filt)+1*FFrac

Assuming you truncate the result.
If FFRac is of another magnitude then change the '1' to
( i.e. an 'halfbit' of the output ).

/Tony

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body  Hi,
just to clarify I added parantheses:

Change:        Filt <-- Filt + (New - Filt)*FFrac

to:         Filt <-- Filt + ((New - Filt)+1)*FFrac

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body   >     some months ago Olin posted the formula that changed DSP world: :)
>         Filt <-- Filt + (New - Filt)*FFrac
>
>     Saddly, my 8 bit implementation have a problem that I can't surpass:
>         filt = filt + ((sig16)new - (sig16)filt) / 2;
>
>     The value never goes to 255 because of integer roundings:
>         filt = 254 + (255 - 254) / 2;
>         filt = 254 + (1)/2;
>         filt = 254;
>
>     There's a way to change this? I've tried to raise the resolution
making:
>         filt = filt + ((sig16)new*2 - (sig16)filt*2 + 1) / 4;
>
>     But got same results because:
>         filt = 254 + ((sig16)255*2 - (sig16)254*2 + 1) / 4;
>         filt = 254 + (510 - 508 + 1) / 4;
>         filt = 254 + (510 - 508 + 1) / 4;
>         filt = 254 + (1) / 4;
>

It helps to have enough fraction bits so that 1/FFrac is not expressed as 0.

*****************************************************************
Embed Inc, embedded system specialists in Littleton Massachusetts
(978) 742-9014, http://www.embedinc.com

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body   Hello Dave,

> I think you correctly recognized that you need to round rather than
> truncate the result, but you only rounded one term instead of the total:
>
> Original:
>       filt = filt + (new - filt) / 2;
>
> Rounded:
>       filt = ((filt*2 + (new*2 - filt*2) / 2) + 1) / 2;
>       filt = ((254*2 + (255*2 - 254*2) / 2) + 1) / 2;
>       filt = ((508 + (510 - 508) / 2) + 1) / 2;
>       filt = ((508 + 2 / 2) + 1) / 2;
>       filt = (509 + 1) / 2;
>       filt = 510 / 2;
>       filt = 255;

you are right! Dumb me...

Saddly, your solution suffers from a problem I haven't seen when I wrote
the original post.

The filter don't go from 1 to 0:

filt = ((filt*2 + (new*2 - filt*2) / 2) + 1) / 2;
filt = ((1*2 + (0*2 - 1*2) / 2) + 1) / 2;
filt = ((2 + (0 - 2) / 2) + 1) / 2;
filt = ((2 + (-1)) + 1) / 2;
filt = (1 + 1) / 2;
filt = 2 / 2;
filt = 1;

Thank you,

Brusque

-----------------------------------------------------------------
Edson Brusque                 C.I.Tronics Lighting Designers Ltda
Research and Development               Blumenau  -  SC  -  Brazil
Say NO to HTML mail                          http://www.citronics.com.br
-----------------------------------------------------------------

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body   Hello Tony,

>Change:     Filt <-- Filt + (New - Filt)*FFrac
>to:         Filt <-- Filt + ((New - Filt)+1)*FFrac
>Assuming you truncate the result.
>If FFRac is of another magnitude then change the '1' to
>( i.e. an 'halfbit' of the output ).

sorry but I don't get it. FFrac have to be between 0 and 1 right? (x/2)
is the same as (x*0.5). So I don't know what you mean with FFrac=4, maybe
(new-filt)/4?

Your formula is great. It overcomes the problem of having to cast New
and Filt variables to 16 bit, but it suffers from the same problem as
Dave's.

Filt <-- Filt + ((New - Filt)+1)/2
Filt <-- 1 + ((0 - 1)+1)/2
Filt <-- 1 + (0)/2
Filt <-- 1

My implementation solved this problem and don't have to deal with
negative numbers:

if (filt < new) filt = filt + ((new - filt)+1) / 2;
else filt = filt - ((filt - new)+1) / 2;

The only problem I've found if when going from 0 to 255 or from 255 to
0:

filt = filt + ((new - filt)+1) / 2
filt = 0 + ((255 - 0)+1) / 2    <--- 255+1=256, it will roll to 0
filt = 0 + (0) / 2
filt = 0 <--- should be 255

Anyway, I can solve this doing 16 bit operations or (easier) with a
couple of "if" clauses.

Thank you very much,

Brusque

-----------------------------------------------------------------
Edson Brusque                 C.I.Tronics Lighting Designers Ltda
Research and Development               Blumenau  -  SC  -  Brazil
Say NO to HTML mail                          http://www.citronics.com.br
-----------------------------------------------------------------

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body   Brusque --

> Saddly, your solution suffers from a problem I haven't seen when I wrote
> the original post.
>
> The filter don't go from 1 to 0:

Yes, that's rounding for you -- it either biases toward 0 (truncation) or
toward +infinity (adding 0.5 and truncating). If you really need to go all
the way from 0 to 255, then you need something more sophisticated.

"Round to even", as used in some floating-point systems, won't work for
you; you'll still get stuck at 254 at the high end.

One approach would be to "dither" the rounding by simply adding 1 or adding
0 on alternate iterations.

-- Dave Tweed

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv mitvma.mit.edu with SET PICList DIGEST in the body  >One approach would be to "dither" the rounding by simply adding 1 or adding
>0 on alternate iterations.

I think that there is no way to make that low pass function span the
entire output domain. Because the math specifies an asymptotical approach
of the output value to the input value and asymptotes do not exist in the
integer math world. I also think that the fiddling with conditional
clauses will introduce points of discontinuity into the mapping function
and if this is a part of a servo or some control loop it may react
strangely to those points.

Peter

--
http://www.piclist.com hint: To leave the PICList
piclist-unsubscribe-request mitvma.mit.edu   > >One approach would be to "dither" the rounding by simply adding 1 or
> >0 on alternate iterations.
>
> I think that there is no way to make that low pass function span the
> entire output domain. Because the math specifies an asymptotical approach
> of the output value to the input value and asymptotes do not exist in the
> integer math world.

True, which is why a few fraction bits do allow the filter to span the
entire range.  I usually make the filter value fixed point with a few extra
bits on the right.  For example, a common case is filtering a 10 bit A/D
value.  I use 16 bits for this with the 10 bit A/D value in the high bits.
That leaves 6 fraction bits for filter roundoff.  It will work exactly as
long as the filter fraction is 2 ** -N, where N is from 0 to 6.

By the way 1/64 is a very heavy filter for most purposes because it has such
a slow response time.  It has a 50% step response time of 45 iterations and
90% in 147 iterations.  Splitting this into two poles of 1/8 each increases
the cycles to compute it but also the response time.  The 50% response drops
to 12 iterations, and the 90% response to 29 iterations.  Both filters have
the same 18dB attenuation of random noise.

> I also think that the fiddling with conditional
> clauses will introduce points of discontinuity into the mapping function
> and if this is a part of a servo or some control loop it may react
> strangely to those points.

I agree.  That kind of fiddling will only lead to trouble.

*****************************************************************
Embed Inc, embedded system specialists in Littleton Massachusetts
(978) 742-9014, http://www.embedinc.com

--
http://www.piclist.com hint: To leave the PICList
piclist-unsubscribe-request mitvma.mit.edu   On Sat, 15 Jun 2002, Olin Lathrop wrote:

> > >One approach would be to "dither" the rounding by simply adding 1 or
> > >0 on alternate iterations.
> >
> > I think that there is no way to make that low pass function span the
> > entire output domain. Because the math specifies an asymptotical approach
> > of the output value to the input value and asymptotes do not exist in the
> > integer math world.
>
> True, which is why a few fraction bits do allow the filter to span the
> entire range.  I usually make the filter value fixed point with a few extra
> bits on the right.  For example, a common case is filtering a 10 bit A/D
> value.  I use 16 bits for this with the 10 bit A/D value in the high bits.
> That leaves 6 fraction bits for filter roundoff.  It will work exactly as
> long as the filter fraction is 2 ** -N, where N is from 0 to 6.

Actually, not entirely true. As I showed in an earlier post, all one needs
to do is reconstruct the equation so that the asymptotes converge in
proper direction. Most of the time when these integer weighted filters are
presented, the new value of the filter is the old value of the filter plus
a term that is weight by an input:

filter = filter + (filter - input)*weight           ....... eq(1)

I agree when constructed this way that it's not possible using integer
arithmetic constrained to the dynamic range of the data to span the entire
dynamic range. However, if the equation is rearranged:

filter = input - (input - filter) * weight'         ....... eq(2)

it is. Mathematically these are identical (weight' is equal to 1-weight).
However, there are several subtle differences when integer arithmetic is
used. In both cases, if the difference the new input and the current state
of filter is less than the reciprocal of the weighting fraction then the
integer arithmetic truncates the result:

if( abs(filter - input) < 1/weight ) then
abs(filter - input) / weight is zero.

In eq(1), this means that the filter will not change states. In eq(2) this
means that the new filter state will become the same as the input. In
other words, if the input spans to the extremes of the dynamic range then
the filter output will too.

Round to Zero
-------------

Now, you don't get something for nothing. As I showed in the PIC assembly
implementation of this algorthm you do need to control the rounding.
Specifically, you need to monitor the sign of the difference and round
towards zero. If the difference is positive, then the integer division
automatically will round toward zero. If the difference is negative you
need to add 1 to the final result. In psuedo code:

; compute   (input - filter) / W
; where input, filter, and W are all integers.
; implement "rounding toward zero"

result = (input - filter) / W

if( result < 0)
result++;

{Quote hidden}

In the general case, this is true. However, the "rounding toward zero"
algorithm has a discontinuity at zero. Fortunately, this is also the point
where the input has no effect on the filter. It's kind of like sin(x)/x
for x equal to zero.

For a weight factor of 1/2, it's interesting to see how the two filter
constructions behave:

(A - B)/2

A-B  round to zero    floor
------------------------------------------
0          0          0
1          0          0
2          1          1
3          1          1
-1          0         -1
-2          0         -1
-3         -1         -2
-4         -1         -2

In other words, the two are the same when the difference is positive. When
the difference is negative, the "round to zero" algorithm is one greater.
The discontinuity is probably easier to see:

rtz      -1 -1  0  0  0  0  1  1  2
floor    -2 -2 -1 -1  0  0  1  1  2
diff     -4 -3 -2 -1  0  1  2  3  4

The round to zero algorithm results in zero for 4 values whereas the floor
algorithm results in zero for only two values. This behavior of the rtz
algorithm is clearly discontinuous. However, this dicontinuity works in
your favor when applied to eq(2).

Let's take an example:

weight = 1/2
filter = 8
input  = 10

eq(1):

1)    filter = 8 + (10 - 8)/2 = 8 + 1 ===> 9
2)           = 9 + (10 - 9)/2 = 9 + 0 ===> 9

eq(2):

1)    filter = 10 - (10 - 8)/2 = 10 - 1 ===> 9
2)           = 10 - (10 - 9)/2 = 10 - 0 ===> 10

Now let's go the other direction.

weight = 1/2
filter = 8
input  = 6

eq(1):

1)    filter = 8 + (6 - 8)/2 = 8 - 1 = 7
2)           = 7 + (6 - 7)/2 = 7 - 1 = 6

eq(2):

1)    filter = 6 - (6 - 8)/2 = 6 - 0 = 6

Here's another interesting test
weight = 1/2
filter = 8
input  = 7  (differ by 1)

eq(1):

1)    filter = 8 + (7 - 8)/2 = 8 - 1 = 7

eq(2):

1)    filter = 7 - (7 - 8)/2 = 7 - 0 = 7

Both filters converge, but the round to zero does so more quickly when the
difference between the input and filter is even.

So as far as non-linearities and discontinuities, I think it's clear to
see that the "frequency response" of both these filters is approximately
the same. In all cases, the number iterations required for the two filters
to reach their final states differ by at most one. The round to zero
implementation however has the benefit of always converging toward the
input value.

Scott

--
http://www.piclist.com hint: To leave the PICList
piclist-unsubscribe-request mitvma.mit.edu   I wrote:
> One approach would be to "dither" the rounding by simply adding 1 or adding
> 0 on alternate iterations.

"Peter L. Peres" <plp ACTCOM.CO.IL> wrote:
> I also think that the fiddling with conditional
> clauses will introduce points of discontinuity into the mapping function
> and if this is a part of a servo or some control loop it may react
> strangely to those points.

Olin Lathrop <olin_piclist EMBEDINC.COM> wrote:
> I agree.  That kind of fiddling will only lead to trouble.

Then you guys simply do not understand the concept of "dither", which
has a long tradition of application in both control systems and audio
applications for *linearizing* the response when dealing with quantization
effects.

A simple "square wave" dither of 1/2 LSB seems perfect for this
application, in terms of low computational cost. Random dither would be
more difficult to compute, and would probably have no net benefit for the
overall application.

Anyway, another approach would be to simply avoid the internal quantization
altogether by maintaining the filter state as a 16-bit variable to begin
with. Just take the high byte as the 8-bit output value when you need it.
Something like this:

filter is a 16-bit variable used to maintain the state of the filter
factor is an 8-bit constant whose value is 256 * FFrac

filter = filter + (new - filter>>8) * factor;

Where "filter>>8" is just C-style notation for "the high byte of filter".
This is an 8-bit subtract, an 8*8->16 multiply (which can be optimized to a
16-bit shift if factor is a power of 2), and a 16-bit addition. The only
cost is one more byte of data storage.

-- Dave Tweed

--
http://www.piclist.com hint: To leave the PICList
piclist-unsubscribe-request mitvma.mit.edu   Hello Scott,

> The filter can be written like so:
>
> Filt' = Filt + (New - Filt)*FFrac
>       = Filt(1-FFrac) + New*FFrac
>
<snip>
> BTW, an efficient implementation of this class of recursive, one-pole
> digital filters can be found:
>
> http://www.dattalo.com/technical/software/pic/twist.asm

this is very interesting. I think now I have something to thought about
in the whole weekend. :)

I'll talk to you later when I got good understanding of your twist.asm.

Best regards,

Brusque

-----------------------------------------------------------------
Edson Brusque                 C.I.Tronics Lighting Designers Ltda
Research and Development               Blumenau  -  SC  -  Brazil
Say NO to HTML mail                          http://www.citronics.com.br
-----------------------------------------------------------------

--
http://www.piclist.com hint: To leave the PICList
piclist-unsubscribe-request mitvma.mit.edu  Hi,
late reply, been out of the loop for a while,

anyway, Edson Brusque wrote:
<snip>
>    The only problem I've found if when going from 0 to 255 or from 255
to 0:
>
>    filt = filt + ((new - filt)+1) / 2
>    filt = 0 + ((255 - 0)+1) / 2    <--- 255+1=256, it will roll to 0
>    filt = 0 + (0) / 2
>    filt = 0 <--- should be 255
>
>    Anyway, I can solve this doing 16 bit operations or (easier) with a
>       couple of "if" clauses.
<snip>
Well remember that you actually have 9 bit's in the working register,
you can also use the carry, in this case, assuming the constants you use
and the flow, would give:

filt = filt + ((new - filt)+1) / 2
filt = 0 + ((255 - 0)+1) / 2    <--- 255+1=256, it will roll to 0 ** and
carry will be set **
filt = 0 + (0) / 2              <---- use RRF, it will roll in the carry
as top bit and produce
filt = 128 <--- should be 255 **nope** should be 128 ( ffrac is = 2 )

/Tony

--
http://www.piclist.com hint: To leave the PICList
piclist-unsubscribe-request mitvma.mit.edu

More... (looser matching)
- Last day of these posts
- In 2002 , 2003 only
- Today
- New search...