Re: [slightly OT] Re: Calculation of filter coefficients given characteristic

From: Eric Scheirer (eds@media.mit.edu)
Date: Tue Nov 23 1999 - 10:49:20 EST


Hi Dana,

Thanks for your great email. I'm going to just think out loud
about how these ideas map into SAOL instrument design for
various requirements. As you said, my interest is not in
real-time use, but I'll try to think about both.

>1. Arbitrary Filter Calculations

>
>It takes a lot of CPU cycles to perform the typical filter design
>calculation. Many transcendental computations are usually needed, such as
>trig, log, etc., not to mention square roots and divides, all of which chew
>up a lot of CPU time. Even at the K-rate, you can be looking at a lot of
>cycles, if you are trying to run something in real time.

I wonder how many applications truly need filter re-design to happen
at the k-rate. At least in the one I had in mind, I'm happy to have
the filters designed at the i-rate, or even at orchestra startup.
Past the issue of how often the filter is changing, there's the issue
of zipper noise that comes from changing parameters on the fly.
How stable are traditional (ie, bilinear transform) methods of
filter design as you change the input parameters? For example,
if I smoothly vary the cutoff frequency of a bandpass designed
with ellip() in Matlab, will the pole/zero locations smoothly
vary?

I suspect this is not a requirement in any applications except
real-time synthesis (maybe some control theory things) and so
it wasn't covered in my filter design class. I suspect that
many traditional methods are pretty bad.

Yet, this is really the method I need, because I was exactly trying
to avoid the pain involved in switching back and forth between
Matlab (to design the filter and spit out the SAOL code) and
SAOL (to do the synthesis). And since I'm not working in real-time
and my filters are static, I don't think the inefficiency matters.

>2. Table Lookup
>
>Throughout history, table lookup of filter coefficients has been by far the
>fastest method for filter coefficient calculation. Pre compute the filter
>coefs, and then store them in a table.
>
>The problem is of course that the table can get big. Also, you need tons
of
>tables to cover all of the cases that you might want.

Yes. In fact, this can be pretty much impossible if you want to allow
score-based control over the filter parameters, since you can't possibly
include *all* filters in the orchestra as precomputed tables.

>3. Table lookup with Interpolation
>
>This idea dates back to before the Beatles, I think. What is new, perhaps,
>is that some filters (at least) interpolate incredibly well.

>
> [...]
>
>This actually is a miracle. If you look a little bit at the algebra, it
>starts to make more sense.

Particularly in the complex plane, it's very straightforward, as
the poles and zeros are being interpolated, I think. You're right
that this is a good trick. And a further generalization of this is
some kind of "basis filter" method, where you have a set of
prototypes and use linear combinations to get the desired
filter as a reconstruction from the set of prototypes, which
acts as a basis for "filter space".

3-D audio people use this trick for data reduction on sets of HRTF
filters, I believe. They get the prototype functions from a Karhunen-
Loeve expansion (principal components analysis) of the set of HRTF
vectors itself. Bill Gardner also showed in his WASPAA paper this
year that you can do nearly as well more efficiently by simply choosing
a good subset of the original HRTF set as your basis set.

>I have not tried other filter designs in a lot of detail, but this trick
>will work with a wide variety of filters. I strongly suspect that elliptic
>filters for example, will work just fine. The algorithm is to pick 1, 2 or
>3 parameters of the ellip() filter that you want to expose to the user as
>parameters, and then generate 2, 4 or 8 sets of filter a[], and b[]
>coefficients to load into a table. Then test to see if the interpolated
>filter frames, when plotted against the original, follow the desired
>constraints.

This would be easy to implement in SAOL. The two-parameter case is
shown for simplicity.

   kopcode interpfilter(table a, table b, ksig p[2], ksig order) {

      // not tested yet

      table a_proto1(data,...);
      table a_proto2(data,...);
      table a_proto3(data,...);
      table a_proto4(data,...); // then the same for b

      tablemap a_proto(a_proto1,a_proto2,a_proto3,a_proto4);
      tablemap b_proto(b_proto1,b_proto2,b_proto3,b_proto4);

      table w1_by_p1(data,...);
      table w2_by_p1(data,...);
      ...
      table w4_by_p2(data,...);

      tablemap w_by_p(w1_by_p1,w1_by_p2,w2_by_p1,...,w4_by_p2);

      ksig i,j,x;

      x = 0; while (x < order) {
        tablewrite(a,x,0); tablewrite(b,x,0);

        i = 0; while (i < 4) {
           j = 0; while (j < 2)

             cur = tableread(a,x);
             w = tableread(w_by_p[i*2+j],p[j]);
             tablewrite(a,x,cur + tableread(a_proto[i],x) * w);

             cur = tableread(b,x);
             w = tableread(w_by_p[i*2+j],p[j]);
             tablewrite(b,x,cur + tableread(b_proto[i].x) * w);
             j = j + 1;
           }
           i = i + 1;
         }
         x = x + 1;
      }
   }

You put the prototype filters in a_protox(...) and b_protox(...) and
the projections of the desired filters by parameter weights in the
wi_by_pj(...) tables, and then this function returns the reconstruction
of the desired filters in the a and b tables. It takes 2*i*j*k MACs
per k-cycle, where i is the order, j is the number of basis filters,
and k is the number of parameters.

To finish this, the next piece of work is to design the prototype
filters in Matlab and calculate the projections of the user parameters
in the prototype basis space.

It might be better to analytically compute the weights from the
user parameters if this is possible, too.

Best,

 -- Eric

+-----------------+
| Eric Scheirer |A-7b5 D7b9|G-7 C7|Cb C-7b5 F7#9|Bb |B-7 E7|
|eds@media.mit.edu| < http://sound.media.mit.edu/~eds >
| 617 253 1750 |A A/G# F#-7 F#-/E|Eb-7b5 D7b5|Db|C7b5 B7b5|Bb|
+-----------------+



This archive was generated by hypermail 2b29 : Mon Jan 28 2002 - 11:46:36 EST