Basics of Fast Discrete Sine Transform

I just got started with DCT/DST but I still fail to understand how the fast DST is supposed to work. The idea behind the FFT is rather apparent and there is very intuitive pseudo-code of it on the Wikipedia article on the Cooley-Tukey algorithm.

The DST-I is defined as $$a_N(k) = \sum_{n=0}^{N-1}x_n \sin\frac{\pi(n+1)(k+1)}{N+1}$$ which can be split into $a_N(k) = b_{N\over 2}(k)+a_{N\over 2}(k)$ and $a_N(N-k-1)=b_{N\over 2}(k)-a_{N\over 2}(k)$, thus exploiting the symmetry in a similar fashion as for the DFT in case of Cooley-Tukey, with $$a_{N\over 2}(k)=\sum_{n=0}^{{N\over 2}-1} x_{2n+1} \sin\frac{\pi(n+1)(k+1)}{N+1},$$ $$b_{N\over 2}(k)=\sum_{n=0}^{{N\over 2}-1} x_{2n} \sin\frac{2\pi(2n+1)(k+1)}{N+1}$$ as explained in e.g. in Britanak, Yip, Rao. Discrete Cosine and Sine Transforms: General Properties, Fast Algorithms and Integer Approximations.

However, I don't get it to describe this as a simple, recursive algorithm as my twiddle factors appear to be wrong. My code (prototyping in python...) looks basically like below, however it yields to the wrong result for any N>2 (yes, N=1 is trivial but apparently it is correct for N=2, so I can't be totally off course?). What am I missing? Or is the problem using this approach that $b_{N\over 2}(k)$ actually is the DST-II, thus this cannot be computed this way at all (i.e. need to compute DST-I of odd parts and DST-II of even parts, recursively?). While it is kind of fun figuring out the solution on my own, any hints are greatly appreciated - this has to be described somewhere this simple, or hasn't it?

def We(N,k):
return math.sin(math.pi*2*(k+1)/(N+1))

def Wo(N,k):
return math.sin(math.pi*(k+1)/(N+1))

def dst_fast(x_in,N):
x = x_in[:]
y = [0]*N
if N == 1:
return x
even = dst_fast(x[::2],N//2)
odd = dst_fast(x[1::2],N//2)
for k in range(N//2):
e = We(N,k) * even[k]
o = Wo(N,k) * odd[k]
y[k] = e + o
y[N-1-k] = e - o
return y

-
A quick check on the correctness of a DST-I code is to make use of the fact that DST-I is nearly self-inverting: applying it twice to your data returns your data multiplied by n/2. –  J. M. Aug 18 '10 at 13:05
After checking my copy of Britanak/Yip/Rao, I suppose your first thought in the error was correct. To quote the book: "...the (n-1)-point DST-I is recursively decomposed into an $\frac{n}{2}$-point DST-II and an $\frac{n}{2}-1$-point DST-I." May I suggest using a real Fourier transform algorithm or a fast Hartley transform for computing the DST-I instead? –  J. M. Aug 18 '10 at 13:14
@J. Mangaldan concerning the correctness, I have a naive implementation (O(N^2)) to compare against, as well as results from dst() from MATLAB. I did not want to use real Fourier transform mainly for two reasons: 1. If I am not mistaken, this would require more than twice as much memory because I need to "mirror" the original input as input to the FFT - or some modifications to the implementation of the FFT itself, both of which is undesired. 2. There are so many fast implementations of DCT, this must be somehow applicable to DST as well. And, after all, I want to know what's wrong. –  user8589 Aug 18 '10 at 14:26
In any event, the book did say each stage of the recursion requires that part of the data be treated with DST-I and the other part be treated with DST-II. As for using FFT for the DST, there's a way to transform your data so you won't need to "mirror" your data to conform to odd symmetry. I believe FFTW has references on how this might be done, but you will have to check. –  J. M. Aug 18 '10 at 15:38
This question might be better for Stack Overflow ... this may be more of a programming question than a conceptual question, and the readers of Stack Overflow tend to be better programmers. –  Peter Shor Aug 19 '10 at 0:26