repak shawahb
it's like killing and eating a real live angel

^
   

rsw@jfet.org


blogroll

        
Fri, 29 Feb 2008

white event

Tonight I was playing around with Octave et al studying the frequency content of sequences generated by a 16-bit LFSR as a function of its feedback polynomial. Using a random selection from Philip Koopman's list of feedback polynomials that produce maximal 16-bit LFSRs, I generated the corresponding sequence and plotted its FFT in Octave.

What I found, effectively, is that more terms in the feedback polynomial whitens the sequence rather substantially. For example, a short feedback polynomial (0x8148) produces the following spectrum:

On the other hand, the spectrum for 0xfff6 looks like this:

This suggests that feedback polynomials with lots of terms are better for dithering sequences. The downside of this, of course, is that each term requires an additional XOR. Fortunately, if you use a Galois LFSR the computation delay is independent of the number of terms.

In case you're curious (more or less identical to what's in the Wikipedia article): 

#include <stdio.h>
#include <stdint.h>

int main(int argc, char **argv)
{
        uint16_t poly;
        uint16_t lfsr = 1;
        uint16_t period = 0;

        if (argc > 1) { poly = (uint16_t) strtoul(argv[1],NULL,0); }
        else { poly = 0xb400; }


        do {
                lfsr = (lfsr >> 1) ^ (-(lfsr & 1) & poly);
                printf("%d\n",lfsr);
        } while (lfsr != 1);
}

All the data, code, et cetera is also available.


[ permalink | 1 comment ]

Mon, 11 Feb 2008

coincident queuage

Over the weekend I finally took the time to code up a little hack to use my AirLinkTek MediaGate MG-35 as the playback mechanism for MP3q. Basically, my solution was to keep the MP3q server running on positron, but to configure the playwrapper.sh script (for those of you who have jacked around with MP3q in the past, you know what I'm talking about) to connect to the MG-35 and start playback.

The default MG-35 firmware has no means by which to do this, but the kind folks over at the mg-35_firmware_mods Yahoo! group have a hacked version that runs a telnet server. Naturally, I'm lazy and don't want to assemble a MIPS uClinux build toolchain, so I just had to make playwrapper.sh smart enough to connect, make sure that positron was mounted, and start playback via telnet. Simple! (A good bit of information about the mp3 player executable on the MG-35 was also gleaned from the mediagate wiki.)

In case any of you want it, you're welcome to the whole MP3q shebang. I'm not really interested in providing tech support for it, so probably if you write to ask questions I'll just berate you for not "getting it." Or maybe just ignore you. Something like that.

As for strange coincidences: my blog is configured such that it always shows the 25 most recent entries on the front page. As a result, each post bumps an old one off the bottom of the list. Strangely, the post that this one is bumping off the post from a bit over a year ago wherein I mentioned that I had found the MG-35 and was intending to make it work with MP3q. A year of good intentions later, it's quite a surprise I'm not already in hell.


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

Sat, 09 Feb 2008

voyage of rediscovery

So it turns out that my discussion yesterday is (as you would expect) a well known one among those who build digital filters. There is even some terminology to go along with it.

The representation of numbers as plus or minus a fractional power of two is effectively canonical signed digit representation, and while my algorithm doesn't explicitly guarantee canonical results, it does so implicitly by minimizing the (log-wise) error in each step. It should be obvious that CSD is more efficient than two's complement at representing numbers, since digits are ternary (0,+,-) rather than binary.


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

Fri, 08 Feb 2008

multiplication is hard, let's go adding

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 Am⋅x[n-m] + Σk=1..L Bk⋅y[n-k]

such that Am and Bk 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 ⋅ 210 ≅ 736 = 0b101101111

Since I'm using fixed point 0.10 notation, the nth 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 ⋅ 210 ≅ 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) ⋅ 2round(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!) ]