Basics of Fast Discrete Sine Transform - MathOverflow most recent 30 from http://mathoverflow.net 2013-05-20T03:45:18Z http://mathoverflow.net/feeds/question/35953 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/35953/basics-of-fast-discrete-sine-transform Basics of Fast Discrete Sine Transform zerm 2010-08-18T11:35:53Z 2011-05-12T05:22:13Z <p>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 <a href="http://en.wikipedia.org/wiki/Cooley-Tukey_FFT_algorithm#Pseudocode" rel="nofollow">Wikipedia article</a> on the Cooley-Tukey algorithm.</p> <p>The DST-I is defined as <code>$$a_N(k) = \sum_{n=0}^{N-1}x_n \sin\frac{\pi(n+1)(k+1)}{N+1}$$</code> which can be split into <code>$a_N(k) = b_{N\over 2}(k)+a_{N\over 2}(k)$</code> and <code>$a_N(N-k-1)=b_{N\over 2}(k)-a_{N\over 2}(k)$</code>, thus exploiting the symmetry in a similar fashion as for the DFT in case of Cooley-Tukey, with <code>$$a_{N\over 2}(k)=\sum_{n=0}^{{N\over 2}-1} x_{2n+1} \sin\frac{\pi(n+1)(k+1)}{N+1},$$</code> <code>$$b_{N\over 2}(k)=\sum_{n=0}^{{N\over 2}-1} x_{2n} \sin\frac{2\pi(2n+1)(k+1)}{N+1}$$</code> as explained in e.g. in <em>Britanak, Yip, Rao. Discrete Cosine and Sine Transforms: General Properties, Fast Algorithms and Integer Approximations</em>.</p> <p>However, I don't get it to describe this as a simple, recursive algorithm as my <em>twiddle</em> 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 <code>$b_{N\over 2}(k)$</code> 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?</p> <pre><code>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 </code></pre> http://mathoverflow.net/questions/35953/basics-of-fast-discrete-sine-transform/35962#35962 Answer by Federico Poloni for Basics of Fast Discrete Sine Transform Federico Poloni 2010-08-18T12:41:10Z 2010-08-18T12:41:10Z <p>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.</p>