I thought I’d share a bit of code that illustrate the expressiveness of
the Csound 7.0 language. With the new additions, it is now possible to
write the code for a radix-2 FFT very compactly using two UDOs
The first one takes care of the reordering of data to do decimation
in time. The second is the in-place transform itself. Both take advantage
of pass-by-reference and the complex data type.
// data reorder for decimation in time (index bit reverse)
opcode Reorder(s:Complex[]):void
N:k = lenarray(s)
i:k, j:k, m:k init 0, 0, 0
while i < N do
if j > i then
tmp:Complex = s[i]
s[i] = s[j]
s[j] = tmp
endif
m = N/2
while m >= 2 && j >= m do
j -= m
m *= 0.5
od
j += m
i += 1
od
endop
// in-place complex-to-complex FFT
opcode FFT(s:Complex[],fwd:b):void
N:k, n:k = lenarray(s), 1
pi:i = fwd == true ? -$M_PI : $M_PI
Reorder(s)
while n < N do
o:k, n2:k = pi/n, n*2
wp:Complex = cos(o), sin(o)
w:Complex = 1,0
m:k = 0
while m < n do
k = m
while k < N do
i:k = k + n
even:Complex = s[k]
odd:Complex = w*s[i]
s[k] = even + odd
s[i] = even - odd
k += n2
od
w *= wp
m += 1
od
n *= 2
od
if fwd then
s = s/N
endif
endop
Here’s an example
n:k, N:k init 0, lenarray(s)
while n < lenarray(s) do
s[n] = cos(2*$M_PI*n/N), 0
n+=1
od
FFT(s,true)
The complex product w *= wp is still in a PR, but the rest is as per the latest beta release.
