Ok, I have made a *little* progress here. Basically I now have a DIF (decimation in frequency) algorithm that is almost complete:

Code:

void dif_transform( int samples[], int length )
{
int i, j, k;
int even_base, odd_base, even, odd;
int block_count, points_per_block, butterflies;
block_count = 1;
points_per_block = 1 << length;
for( i = 0; i < length ; ++i )
{
butterflies = points_per_block >> 1;
even_base = 0;
for( j = 0 ; j < block_count; ++j )
{
odd_base = even_base + butterflies;
for( k = 0; k < butterflies; ++k )
{
even =
samples[ even_base + k ] + samples[ odd_base + k ];
odd =
( samples[ even_base + k ]
- samples[ odd_base + k ] )
* twiddle( points_per_block, k );
samples[ even_base + k ] = even;
samples[ odd_base + k ] = odd;
}
even_base += points_per_block;
}
block_count <<= 1; // twice as many blocks
points_per_block >>= 1; // half as many points in each block
}
return;
}

I don't take any credit for it, it was taken from the example:

Code:

PROCEDURE DIF(p,VAR f);
LOCAL Bp,Np,Np',P,b,n,BaseE,BaseO,e,o;
BEGIN {DIF}
{initialise pass parameters}
Bp:=1; {No. of blocks}
Np:=1<<p; {No. of points in each block}
{perform p passes}
FOR P:=0 TO (p-1) DO
BEGIN {pass loop}
Np':=Np>>1; {No. of butterflies}
BaseE:=0; {Reset even base index}
FOR b:=0 TO (Bp-1) DO
BEGIN {block loop}
BaseO:=BaseE+Np'; {calc odd base index}
FOR n:=0 TO (Np'-1) DO
BEGIN {butterfly loop}
e:= f[BaseE+n]+f[BaseO+n];
o:=(f[BaseE+n]-f[BaseO+n])*T(Np,n);
f[BaseE+n]:=e;
f[BaseO+n]:=o;
END; {butterfly loop}
BaseE:=BaseE+Np; {start of next block}
END; {block loop}
{calc parameters for next pass}
Bp:=Bp<<1; {twice as many blocks}
Np:=Np>>1; {half as many points in each block}
END; {pass loop}
END; {DIF}

My primary problem now is with the "twiddle-factor" function. How to implement it? It was mentioned that a table lookup would be most efficient. Does anyone know about this? Also, since the output is bit-reversed, is there an easy way to reorder the bits efficiently?