[Math] Basics of Fast Discrete Sine Transform

algorithmsfourier analysis

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

Best Answer

Not an expert with DST, but I think you're missing a costant factor when using dst_fast(something,N//2) to build dst_fast(something,N) --- the former contains a N/2+1 denominator, the latter needs a N+1 denominator. Therefore there should be a conversion factor with a term looking like N/2+1, which I don't see in your code.

Related Question