repak shawahb

dry humor for cactus farmers

200812 ( 1)
200811 ( 2)
200809 ( 1)
200808 ( 1)
200806 ( 2)
200805 ( 2)
200804 ( 2)
200803 ( 5)
200802 ( 4)
200801 ( 1)
200712 ( 5)
200711 ( 3)
200709 ( 1)
200706 ( 1)
200704 ( 4)
200703 ( 1)
200702 ( 3)
200701 ( 5)
200612 ( 3)
200610 ( 2)
200609 ( 2)
200608 ( 3)
200607 (10)
200606 (10)
200605 ( 9)
200604 ( 7)
200603 ( 9)
200602 (11)
200601 (26)
200512 (20)
200511 (22)
200510 (19)
200509 (13)
200508 ( 2)
200507 ( 3)
200506 ( 7)
200505 (26)
200504 ( 9)
200503 (15)
200502 (20)
200501 (28)
200412 (16)
200411 ( 5)
200410 (18)
200409 ( 7)
all

blogroll

Fri, 08 Feb 2008

I was considering this problem the other day, and while it's simple to arrive at a solution, I rather like my particular solution. Thus:

Let's say that you want to implement a digital filter of the form

y[n] = Σ_m=0..N A_m⋅x[n-m] + Σ_k=1..L B_k⋅y[n-k]

such that A_m and B_k are arbitrary precision constants. An example is this 2nd order Chebyshev(1) filter:

y[n] = .0027417⋅x[n] + .0054834⋅x[n-1] + .0027417⋅x[n-2] + 1.90020⋅y[n-1] - .91970⋅y[n-2]

(this filter is the result of the octave cheby1() function). An obvious way of doing this is to build some number of multipliers and adders to compute each term. But it's also possible to do this with only an adder as long as we are allowed to take time to do it. The way we do this is by expressing the constant coefficients as sums of powers of two and then do shifts and addition on the variables. For example,

.75⋅x = x - .25⋅x = x - x>>2

.75 is obviously a contrived example, but in general how do I come up with a representation of a number as such a sum? The simplest way I can think of is to express the number in fixed point binary with a chosen precision and then use the bits as a guide. Consider computing .718282 (the fractional portion of e) to 10 bits.

.718282 ⋅ 2¹⁰ ≅ 736 = 0b101101111

Since I'm using fixed point 0.10 notation, the n^th bit from the left represents 2^-n. The bits that are 0 are effectively the missing fractional pieces, so I quickly surmise that

.718282 ≅ (1 - 1/1024) - 1/4 - 1/32

(The reason I start with (1 - 1/1024) is because 1 is not actually representable in 0.10 notation; the largest possible number is 1 - 1/1024. If you ignore this subtraction, you will be accurate to 9 bits instead of 10, so you can always just produce the result you want by computing to one more bit than necessary and dropping this correction factor.)

This works great! Let's do another one, say, π/4.

π/4 ≅ .7853982 .7853982 ⋅ 2¹⁰ ≅ 804 = 0b1100100100 ∴ π/4 ≅ (1 - 1/1024) - 1/8 - 1/16 - 1/64 - 1/128 - 1/512 - 1/1024 = 1 - 1/8 - 1/16 - 1/64 - 1/128 - 1/256

Ugh! That is a helluva lot of subtractions. Why so many? Note that .7853982 is just slightly greater than .75, and thus a better way to approximate this would have been to subtract .25 and then start adding small correction factors instead of sticking strictly to subtraction. What we've uncovered is that our initial algorithm, while yielding correct results, is very suboptimal.

How can we improve it? Well, in order to get to the desired accuracy in the fewest number of steps, clearly we want to take the biggest chunk possible out of the error term. This implies that we want an iterative algorithm that makes a guess, then computes the error, then approximates the error as a power of two and recomputes the approximation:

do { error = estimate - number eapprox = sign(error) ⋅ 2^{round(log2(abs(error)))} estimate = estimate - eapprox } while ( error > tolerance )

Now, running this algorithm to the same precision on π/4 yields

π/4 ≅ 1 - 1/4 + 1/32 + 1/256

I'm sure there's some way of doing this from the fixed point trickery above as well. Let's call that an exercise for the reader.

[ permalink | 0 comments (add one you lazy bastard!) ]

writebacks (add one you lazy bastard!)

post a comment: