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
```