astomo.us

Alvin Lucier’s “I am Sitting in a Room” is a classic piece of feedback-based sound design. It consists of Lucier sitting in a room with 2 distant mics, 2 tape decks, and speakers. The recording begins with Lucier verbally describing the process of the piece, which is that he is recording his voice, and after recording, he will play that back through the speakers, and record that through the mics onto the 2nd tape deck. This continues for multiple generations. At first, you just hear reverberation of the rooms acoustics as he replays the tape into the room. Gradually, the resonant frequencies of the room (i.e. sinusoids whose wavelength divides the distance from speaker to mic) begin to stand out, and by the middle of the 2nd side of the lp the original recording is completely obscured.

Last month I wrote a max patch that implements a direct echo/reverb feedback loop to digitally simulate Lucier’s design, with a few extra features. There are several presets you can mess around with- the latter preset object implements a simulated acoustic reverb network from James Anderson Moorer’s 1978 paper ‘about this reverberation business’. Ive also included an example aiff file and its output after several passes through the loop.

NOTE OF CAUTION- if you start toying with the reverb network- the cascade output gain isnt normalized so the feedback tends to get out of control after a few iterations if you dont manually lower the gain so mind your ears.

Perhaps you werent aware that there is a GIANT HEXAGON ON SATURN’S NORTH POLE. It appears to be a variety of standing wave. These are new pictures with an unprecedented level of detail- note the shockwaves (at least ithink theyre shockwaves) propagating off of the corners. Thanks to Joe Mason from the Cassini imaging media team for fwding these to me.

Two thoughts, in particular, stand out. In any sufficiently large field of study:

(a) there is no single person who is an authority, and

(b) at any given moment, there will be many people working independently on the same problem. To anyone who has done any type of research, these are self-evident truths. Ignorance and uncertainty, as well as a vast and chaotic literature are inevitable by-products of a free research environment. When there is no such freedom, it is easy for one to become an authority and, once an authority, to make sure that everyone in one’s charge is working together on the same problem; at that stage, it is only a small step to further ensure, as Lysenko did in Soviet Russia, that everyone uncovers the same results, their truth notwithstanding.

When presented with an Aphex Twin track, Karlheinz Stockhausen once commented: “I heard the piece Aphex Twin of Richard James carefully: I think it would be very helpful if he listens to my work “Song of the Youth,” which is electronic music, and a young boy’s voice singing with himself. Because he would then immediately stop with all these post-African repetitions, and he would look for changing tempi and changing rhythms, and he would not allow to repeat any rhythm if it [was] varied to some extent and if it did not have a direction in its sequence of variations.”

To which AT responded: “I thought he should listen to a couple of tracks of mine: ‘Didgeridoo,’ then he’d stop making abstract, random patterns you can’t dance to”.

- “Advice to Clever Children”. The Wire 67 (141): 553.

The other day I was listening to AT’s “Selected Ambient Works Vol. 2″ and track #2 off the first disc got me thinking. In case you’ve never heard of R.D. James (aka Aphex Twin), he is an electronic musician. Occasionally he does something nerdy and/or snarky which cracks me up. For example his 1999 single Windowlicker had several images embedded as sonograms inside the audio. Much of his music sounds meticulously edited, but occasionally there are melodies or modulation parameters or something that sound algorithmic. For fun I decided to try and ‘cover’ track #2 algorithmically with no editing. To construct the melody and modulation data I used a bit-reversal algorithm on an ascending harmonic scale.

Bit-reversal (pictured above as appiled to an 8-element array) is, as the name suggests, an operation (permutation actually) wherein you take an array of data and map each element to the unique address whose binary digits are the mirror image of those in its current address. It is a low-level primitive that comes up surprisingly often in computational harmonic analysis. The most prominent example is the FFT of course, but there are several other interesting transforms (such as the sequency ordered Walsh-Hadamard transform) that make similar use of it. Bit-reversal is plenty interesting on its own however. It can be done in-place and in O(n) time and can be profitably parallelized on SIMD architectures (such as GPU’s). To create my data I used a recursive version (written in Python):

from numpy import *

def bit_reverse_traverse(a):
    n = a.shape[0] 
    if n == 1:
        yield a[0]
    else:
        even_index = arange(n/2)*2
        odd_index = arange(n/2)*2 + 1
        for even in bit_reverse_traverse(a[even_index]):
            yield even
        for odd in bit_reverse_traverse(a[odd_index]):
            yield odd
    return

def get_bit_reversed_list(l):
    n = len(l)
    index = arange(n)
    b = []
    for i in bit_reverse_traverse(index):
        b.append(l[i])      
    return b

The first function bit_reverse_traverse is a generator, which is a sortof functor (by which I mean function with internal state, not the category-theoretic kind) that behaves like an iterator. However instead of creating the entire array in memory and feeding it out one element at a time, generators rely on internal state (using the yield statement, which acts like a stack frame freezer) to synthesize a new element each time their next method is called. Hence they are much more efficient than normal iterators (at the cost of not being random access) and enable one to do things like bit-reverse a 65,536 element array on a laptop reasonably quickly. This generator works a bit like merge-sort: by recursively iterating over copies of itself applied to 2 sub-arrays of even and odd-numbered elements (accessible using numpy syntax a[even_index]) until the arrays are of size 1, then it just returns the element. The second function get_bit_reversed_list actually performs the bit-reversal on a list by iterating over the list using the generator and appending each array element returned by the yield statements.

One can do all sorts of interesting things with generators- usually with some kindof explicit recursion (indeed they are conceptually quite similar to basic streams in Scheme), but not always. For example you can easily phrase the Sieve of Erastosthenes algorithm for generating prime numbers as a generator:

from numpy import any

def primes():
    n = 2
    p = []
    while True:
        if not any( n % f == 0 for f in p ): 
            yield n
            p.append( n )
        n += 1

As you can see the code isn’t recursive.

Anyways, after creating some lists this way I converted them into MIDI data and played them on a sampling synthesizer. Here is the result.

Today I was reading the slides from Felix von Leitner’s talk on C compiler optimizations (given at the recent Linux Kongress). Basically his point is that there is no point in making your code confusing in an attempt to optimize it b/c usually the resulting object code will not be faster. I’m not sure how true that is in general (it doesn’t seem to be the case with OpenCL at least in my experience), but I really like his point. The slides are also a good read if you like learning semi-obscure C tricks. For example if you want to inline a function in C its best to actually not use the inline keyword. This ‘destroys code locality’ according to Felix which Im assuming means that it can cause pointers to what is apparently the same function to compare not equal to one another. However what you can do is declare them static if you’re calling them exactly once. In that case the compiler will inline the function. For example, when compiled with the command gcc -S -O2, the code:

static long abs2(long x) {
return x>=0?x:-x;
} /* Note: > vs >= */

long abs3(long x) {
return x>=0 ? x : -x;
}

int main(int argc, char *argv[]) {
int x = 11;
int y;

y = abs2(x);
y = abs3(x);

return 0;

}

will produce assembly that contains no definition of abs2. However if you add a second call to abs2 and look at the assembly again you’ll find that all of a sudden you have an extra:

_abs2:
LFB2:
  pushq  %rbp
LCFI0:
  movq  %rsp, %rbp
LCFI1:
  movq  %rdi, %rdx
  sarq  $63, %rdx
  movq  %rdx, %rax
  xorq  %rdi, %rax
  subq  %rdx, %rax
  leave
  ret
LFE2:
  .align 4,0x90

When he was told that Apple Computer had just bought a Cray to help design the next Apple Macintosh, Cray commented that he had just bought a Macintosh to design the next Cray.

- According to Jim Gray (quoted by C. Gordon Bell in his “Seymour Cray Perspective”)

SIAM is holding a conference on parallel computing in Seattle next February. I’m pretty psyched about it mostly because I’ve been sinking a lot of time into GPU-based stream computing lately and from a cursory look at the program I’d say that GPGPU is going to be the subject of a lot of interest.

For example here is an abstract from one of the Joint JSIAM — SIAM Lectures (sortof amusingly entitled: GPU Acceleration: A Fad or the Yellow Brick Road onto Exascale?):

“Since the first commodity x86 cluster Wigraf achieving paltry 10s~100s Megaflops in 1994, we have experienced several orders of magnitude boost in performance. However, the first Petaflop was achieved with the LANL RoadRunner, a Cell-based “accelerated” cluster, and in 2010 we may see the first (GP)GPU-based cluster reaching Petaflops. Do such non-x86“accelerator” merely push the flops superficially, or are they fundamental to scaling? Based on experiences from TSUBAME, the first GPU-accelerated cluster on the Top500, we show that GPUs not only achieve higher performance but also better scaling, and in fact their true nature as multithreaded massively-parallel vector processor would be fundamental for Exascale. Such results are being reflected onto the design of TSUBAME2.0 and its successors.”

In general the computer hardware industry is unpredictable and early adopters often get burned. How widespread GPGPU eventually becomes still remains to be seen, but such talk must be encouraging to anyone who has invested time into constructing and implementing algorithms on GPU architectures. Anyways, I’m hoping to make the trip up and do a poster presentation of some of the GPU-based matrix compression work I’ve been doing.

In my last post I made an off-hand remark about writing tail-recursive code in Python. I thought that I would elaborate on this. Python’s creator, Guido van Rossum, has famously said that he’s not interested in implementing tail recursion in the language- much to the dismay/disgust of functional programmers everywhere. Im certainly not a hard-core functional programmer and hadnt devoted much thought to this until a few weeks ago when my friend Judah gave a talk on functional programming in numerical analysis; he made a pretty convincing case that using a functional style has a lot of advantages. Of course there’s not much point in doing this if your language doesnt do tail-recursion optimization b/c your code will most likely explode the stack and/or run abysmally slowly. Im still not sure about the second problem, but one (admittedly hacky) way to avoid the first problem when running tail-recursive code in Python is to use decorators (not to be confused with the Decorator design pattern).

So what are decorators in Python? Well, one way to think of them is as a Pythonic analog to the classic C preprocessor macro: they modify existing code at (byte) compile time. However decorators have a somewhat more functional flavor, for example:

#decorators.py
import sys

class simple_decorator(object):
  def __init__(self, f):
    print "initializing simple_decorator"
    f() # Prove that function definition has completed

  def __call__(self):
    print "calling simple_decorator"

def foo1():
  print "calling foo1"

@simple_decorator
def foo2():
  print "calling foo2"

foo1()
foo2()

You can see that simple_decorator is a class whose constructor takes a function as its second argument (functions are first-class objects in Python). The call method typically changes the behavior of the base function in some way. Here it completely replaces it. So in my python interpreter:

In [5]: from decorators import *
initializing simple_decorator
calling foo2
calling foo1
calling simple_decorator

In [6]: foo1()
calling foo1

In [7]: foo2()
calling simple_decorator

foo1 and foo2 have the same source code (modulo the print string of course), but different behavior- thats because foo2 was modified using simple_decorator. We can use this construction to implement some really basic tail-recursion ‘optimization’ by collecting local vars each time the stack explodes and piping them back into the function:

class stack_explosion:
  def __init__(self, args, kwargs):
    self.args = args
    self.kwargs = kwargs

def tail_recursive(foo):
  def foo_tr(*args, **kwargs):
    f = sys._getframe()
    if f.f_back and f.f_back.f_back and f.f_back.f_back.f_code == f.f_code:
      raise stack_explosion(args, kwargs)
    else:
      while True:
        try: return foo(*args, **kwargs)
        except stack_explosion, boom:
          args = boom.args
          kwargs = boom.kwargs
  return foo_tr

@tail_recursive
def fact1(n, result=1):
  if n == 0:
    return result
  return fact1(n-1, n*result)
  
def fact2(n, result=1):
  if n == 0:
    return result
  return fact2(n-1, n*result)

print fact1(10000)
print fact2(10000)

Running fact1(10000) returns 10000! as one would hope, wheras fact2(10000) (the undecorated version) returns a runtime stack-overflow error. Cool, right? This approach was originally pointed out to me by my friend John. Actually speeding up tail-recursive code in Python is more complicated though, and depends a lot on which implementation youre using. One way to go, at least in CPython, is to use a trampoline. This (I believe) is how Scheme does it.

Quick edit: According to Craig you use trampolines to implement tail recursion in scheme if you’re building it on top of a system that doesn’t handle tail calls correctly. If you’re just compiling to assembly, though, you don’t need a trampoline: CPS the program and you can correctly generate the jumps when you’re emitting the assembly.

Ive just finished reading Godfried Toussaint’s excellent paper ‘The Euclidean Algorithm Generates Traditional Musical Rhythms’.   As the title suggests, Toussaint explores the applications of Euclid’s famous algorithm for computing the GCD of two integers to algorithmic rhythm generation.  It turns out that if you take a templated version of the algorithm and apply it to list types in place of integers with a certain mix operation in place of subtraction, and you feed it lists of the type [1,1,1.....,0,0....,0],  you get a procedure that tries to distribute the zeros amongst the ones with maximal uniformity.  For example, applying the procedure to the list [1,1,1,1,0,0,0,0] returns [1, 0, 1, 0, 1, 0, 1, 0] as you might expect.  Things get more interesting when you introduce some non-uniformity however.  Applying it to [1,1,1,0,0,0,0,0] returns [1, 0, 0, 1, 0, 1, 0, 0], a common clave known as tresillo.

The procedure I’ve just described was first applied to pulse timing in spallation neutron source (SNS) accelerators by Eric Bjorklund, but the basic premise of trying to spread stuff out as much as possible subject to certain constraints pops up everywhere in mathematics, from the distribution of zeros of the zeta function to the eigenvalues of a Gaussian ensemble (the two are actually related, a fact famously discovered in the cafeteria of the Institute for Advanced Study) to Ornstein-Uhlenbeck processes (O-U processes, incidentally, are the continuous-time version of what you get when you feed pure white noise into a certain kind of IIR filter).   Not to mention various elementary physics phenomena such as the diffusion equation and magnetism.  This stuff is still an active area of research in mathematics.  If youre interested I suggest checking out Terry Tao’s blog on the subject for starters.

Anyways, this is incredible and I could go on and on about it, but its much more fun to just implement Bjorklund and start messing around.  I find scripting languages in general and Python in particular to be a good tool for this sort of thing. The remainder of the post is a sortof walk-through of my python code, so fire up your Python interpreter and read on.

For reference, here is a subtractive version of the EA:

def gcd(a,b):
  if a==0: return b
  while b!=0:
    if a>b: a = a-b
    else: b = b-a
  return a

To morph EA land into Bjorklund we need a analogs of the integers a and b, >, and -. In Bjorklund a and b are lists and are at times glued together into one list. So we need a way of splitting a list into two components a and b. The split works by separating out the largest homogeneous chuck of the tail of the list. I did this in a tail recursive style:

def split(a,b=[]):
  b.insert(0,a.pop())
  if not a or b[0] != a[-1]: return a,b
  else: return split(a,b)

So for example in my python interpreter:

In [64]: a = [[1],[1],[1],[0],[0],[0],[0],[0]]

In [65]: split(a)
Out[65]: ([[1], [1], [1]], [[0], [0], [0], [0], [0]])

In [66]: a = [[1, 0], [1, 0], [1, 0], [1, 0]]

In [67]: split(a)
Out[67]: ([], [[1, 0], [1, 0], [1, 0], [1, 0]])

In the latter case the entire list is homogeneous so split returns an empty list and the original list. This will be our base case (ie the Bjorklund version of b==0).

The analog of subtraction is mix.

def mix(a,b):
  m = min(len(a),len(b))
  if len(a) > m:
    return [a[i]+b[i] for i in range(m)] + a[m:]
  else:
    return [a[i]+b[i] for i in range(m)] + b[m:]

Mix takes two lists and concatenates them element-wise, then concatenates the result with whatever is left over. For example:

In [72]: a = [[1],[1],[1],[0],[0],[0],[0],[0]]

In [73]: b = [[1],[1],[1],[1],[1],[0],[0]]

In [74]: mix(a,b)
Out[74]: [[1, 1], [1, 1], [1, 1], [0, 1], [0, 1], [0, 0], [0, 0], [0]]

Now the morphed EA looks like this:

def foo(a):
  x,y = split(a)
  while len(x) != 1 and len(x) != 0:
    x,y = split(mix(x,y))
  return [x,y]

One final primitive we need is a function to flatten out the result:

def flatten(x):
  result = []
  for el in x:
    if hasattr(el, "__iter__") and not isinstance(el, basestring):
      result.extend(flatten(el))
    else:
      result.append(el)
  return result

Finally, the Bjorklund algorithm can be written succintly as:

def bjorklund(a):
  return flatten(foo(a))

So for example to create a tresillo clave and its ‘mirror image’:

In [78]: a = [[1],[1],[1],[0],[0],[0],[0],[0]]

In [79]: bjorklund(a)
Out[79]: [1, 0, 0, 1, 0, 1, 0, 0]

In [80]: b = [[1],[1],[1],[1],[1],[0],[0],[0]]

In [81]: bjorklund(b)
Out[81]: [1, 0, 1, 1, 0, 1, 0, 1]

Of course this is just the tip of the iceberg. As a white rabbit, observe that the operator is contravariant (i.e. bjorklund(a+b) = bjorklund(b) + bjorklund(a) where + here denotes list concatenation):

In [82]: a = [[1],[1],[1],[0],[0],[0],[0],[0]]

In [83]: b = [[1],[1],[1],[1],[1],[0],[0],[0]]

In [84]: bjorklund(a+b)
Out[84]: [1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1]

Addendum: If Python isnt your thing, there are a few other implementations of Bjorklund out there in Ruby and Lisp.