I'm trying to implement Algorithm 2 from this paper:

Xuancheng Shao and Steven G. Johnson: Type-IV DCT, DST, and MDCT algorithms with reduced numbers of arithmetic operations, https://arxiv.org/pdf/0708.4399.pdf

It outlines an algorithm for computing an even-sized DCT type 4 via a DCT type 3 and a DST type 3 of half size.

I've made a few attempts at implementing this, but either the paper has a mistake or I have a serious misunderstanding about what's happening here.

https://pastebin.com/uFiHsMKX Here's my C++ code, which uses FFTW to compute the inner transforms, as well as to compute an "expected output". Aside from FFTW all the code necessary to compile is in this one file. The function compute_dct4 is what does the actual processing -- I did my best to name variables the same as in the paper, so it should be clear what in the code corresponds to what in the paper.

I'd appreciate some extra eyes on what's going wrong. Thanks!

  • $\begingroup$ I hope you are aware of the fact that the proposed algorithm is at most 5% more efficient than the previous? published algorithms. And also why do you need DCT-IV and not DCT-II or DCT-I which are the most common types for multimedia applications ? $\endgroup$
    – Fat32
    Apr 3, 2017 at 23:05
  • $\begingroup$ I'm writing a general-purpose DCT library. This library will have DCT and DST types 1-4. It will also make the MDCT available, which is built on top of the DCT-4 $\endgroup$
    – Elliott M
    Apr 4, 2017 at 2:23
  • $\begingroup$ I'm not trying to implement their extra-special algorithm, I just used this paper as a reference because their description of the original algorithm is concise and readable. I've read the paper that originally describes the algorithm, and it doesn't provide any insight into what's going wrong, sadly. $\endgroup$
    – Elliott M
    Apr 4, 2017 at 2:23

1 Answer 1


I found the solution. There was a mistake in the paper.

On the first 2 lines of Algorithm 2:

w0 ← x0
vN/2 ← xN−1

x0 and xN-1 should be scaled by 2, giving:

w0 ← x0 * 2
vN/2 ← xN−1 * 2

The following paste has working code (and an open source license, so feel free to copy it)



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.