Derivation of fast DCT-4 algorithm based on DFT

2 Replies

It’s well known that an N point DCT-4 can be computed using an N/2 point complex FFT. Although the algorithm is widespread, the texts which I have read on the subject have not provided the details as to how it works. I’ve been trying to understand this for a while (I tend not to use algorithms until I understand them) and finally figured it out and thought I would share (please provide comments/links if you can find a better or shorter explanation – this is as good as I could get).

Start with the definition:

    \[ X_k = \displaystyle\sum\limits_{n=0}^{N-1} x_n \cos \frac{ \pi \left( n + \frac{1}{2} \right) \left( k + \frac{1}{2} \right) }{N} \]

Trig functions are hard to manipulate in expressions (for me anyway), so knowing that x_n is real, we change the \cos into a complex exponential and obtain the following:

    \[ X_k = \Re \left\{ \displaystyle\sum\limits_{n=0}^{N-1} x_n \mathrm{e}^{\frac{\mathrm{j} \pi \left( n + \frac{1}{2} \right) \left( k + \frac{1}{2} \right) }{N}} \right\} \]

It is worth noting that the sign of the exponent is irrelevant.

We then split the expression up to operate over two half-length sequences composed of x_{2n} and x_{N-1-2n}. The second sequence is reversed and decimated because when the n is replaced by N-1-2n in the n + \frac{1}{2} term of the exponential, the expression is negated rather than the offset being modified. We move the N term into another exponential which becomes trivial as N cancels the denominator. i.e.:

    \[ X_k = \Re \left\{ \displaystyle\sum\limits_{n=0}^{N/2-1} x_{2n} \mathrm{e}^{\frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( k + \frac{1}{2} \right) }{N}} + x_{N-1-2n} \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( k + \frac{1}{2} \right) }{N}} \mathrm{e}^{\mathrm{j} \pi \left( k + \frac{1}{2} \right)} \right\} \]

Or:

    \[ X_k = \Re \left\{ \displaystyle\sum\limits_{n=0}^{N/2-1} x_{2n} \mathrm{e}^{\frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( k + \frac{1}{2} \right) }{N}} + \mathrm{j} {\left( -1 \right)}^{k} x_{N-1-2n} \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( k + \frac{1}{2} \right) }{N}} \right\} \]

The exponentials still contain no term which resembles an DFT (over length N/2 anyway). To get one, we break up X_k into the terms X_{2k} and X_{N-1-2k}. Again we choose to reverse the second sequence because it will keep the exponential terms in a similar but negated form.

    \[ X_{2k} = \Re \left\{ \displaystyle\sum\limits_{n=0}^{N/2-1} x_{2n} \mathrm{e}^{\frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( 2k + \frac{1}{2} \right) }{N}} + \mathrm{j} x_{N-1-2n} \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( 2k + \frac{1}{2} \right) }{N}} \right\} \]

    \[ X_{N-1-2k} = \Re \left\{ \displaystyle\sum\limits_{n=0}^{N/2-1} \mathrm{j} x_{2n} \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( 2k + \frac{1}{2} \right) }{N}} - x_{N-1-2n} \mathrm{e}^{\frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( 2k + \frac{1}{2} \right) }{N}} \right\} \]

Now because we are only interested in the real terms, we can ignore or conjugate anything which contributes to the imaginary term of the above expressions. This means that we can conjugate the exponentials when they are only modulating a real term. Doing this, we obtain:

    \[ X_{2k} = \Re \left\{ \displaystyle\sum\limits_{n=0}^{N/2-1} \left( x_{2n} + \mathrm{j} x_{N-1-2n} \right) \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( 2k + \frac{1}{2} \right) }{N}} \right\} \]

    \[ X_{N-1-2k} = \Re \left\{ \displaystyle\sum\limits_{n=0}^{N/2-1} \left( \mathrm{j} x_{2n} - x_{N-1-2n} \right) \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( 2k + \frac{1}{2} \right) }{N}} \right\} \]

Almost there, but to obtain the full benefit we need the inner terms to be the same. If we now multiply the X_{N-1-2k} term by - \mathrm{j}, we get:

    \[ X_{N-1-2k} = - \Im \left\{ \displaystyle\sum\limits_{n=0}^{N/2-1} \left( x_{2n} + \mathrm{j} x_{N-1-2n} \right) \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2n + \frac{1}{2} \right) \left( 2k + \frac{1}{2} \right) }{N}} \right\} \]

Done. If we expand out the X_{2k} and X_{N-1-2k} terms completely we get:

    \[ X_{2k} = \Re \left\{ \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2k + \frac{1}{2} \right) }{2N}} \displaystyle\sum\limits_{n=0}^{N/2-1} \left( x_{2n} + \mathrm{j} x_{N-1-2n} \right) \mathrm{e}^{- \frac{\mathrm{j} \pi n}{N}} \mathrm{e}^{- \frac{\mathrm{j} 2 \pi n k }{N/2}} \right\} \]

    \[ X_{N-1-2k} = - \Im \left\{ \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2k + \frac{1}{2} \right) }{2N}} \displaystyle\sum\limits_{n=0}^{N/2-1} \left( x_{2n} + \mathrm{j} x_{N-1-2n} \right) \mathrm{e}^{- \frac{\mathrm{j} \pi n}{N}} \mathrm{e}^{- \frac{\mathrm{j} 2 \pi n k }{N/2}} \right\} \]

Which give us the steps for a DFT based DCT-4 algorithm:

  • Transform the N point real sequence x_n into the N/2 point complex sequence y_n = x_{2n} + \mathrm{j} x_{N-1-2n}.
  • Multiply each element of the sequence y_n by \mathrm{e}^{- \frac{\mathrm{j} \pi n}{N}}.
  • Find Y, the DFT of the sequence y.
  • Multiply each element of Y_k by \mathrm{e}^{- \frac{\mathrm{j} \pi \left( 2k + \frac{1}{2} \right) }{2N}}.
  • The DCT-4 outputs are given by:

        \[ X_k = \left\{ \begin{array}{l l} \Re \left\{ Y_{k/2} \right\} & \quad \text{for even $k$} \\ - \Im \left\{ Y_{(N-1-k)/2} \right\} & \quad \text{for odd $k$} \end{array} \right. \]

2 thoughts on “Derivation of fast DCT-4 algorithm based on DFT”

  1. Niels Möller

    Note that it’s an easy change to use identical pre- and postfactors,

    e^{- j \pi (k+1/8) / N}

    The interesting property is that when we multiply together all
    the exponentials, the FFT exponent 4nk/N is replaced by

    4 n k/N + (k+1/8) / N + (n+1/8) / N
    = (4n+1)(4k+1) / 4N

    which resembles the DCT “kernel” expression, for indices 2n and 2k.

    Reply ↓

Leave a Reply

Your email address will not be published. Required fields are marked *