Tag Archives: DFT

Three Fast Methods to Compute the DCT-2

After doing my "proof thing" of the fast DFT based implementation of the DCT-4 (which does actually come in useful in my line of work), I thought I would have a fiddle with some of the other DCT variants. The DCT-2 is a transform which I doubt I will ever use, but I managed to find three ways (excluding the slow and obvious one) to build it. Here is the definition:

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

Like in the DCT-4 post, I hate trig functions so the first thing to do is transform  \cos into a complex exponential. I'll use this as the starting point for the methods.

 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) k }{N} } \right\}

Method 1: Recursion and the DCT-4

All of the basis functions of the DCT-2 have a symmetry about  n=\frac{N}{2} ; the symmetry is even for even  k and odd for odd  k . This is a big hint that we can benefit from splitting the input sequence into two half-length sequences where one sequence is reversed:

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

As on the DCT-4 page, we can conjugate any terms we like in the summation as we are only interested in the real part of the result. This allows us to rewrite the above as:

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

Split  X_k into even and odd terms:

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

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

The end. It can be seen that the  X_{2k} terms can be obtained by recursively calling the function again. The  X_{2k+1} terms can be obtained via a DCT-4. Writing out the steps explicitly:

  • Create a sequence  y_n = x_n + x_{N-n-1} of length  \frac{N}{2}
  • Create a sequence  z_n = x_n - x_{N-n-1} of length  \frac{N}{2}
  • The resulting  X_{2k+1} terms are given by the DCT-4 of the sequence  z_n
  • The resulting  X_{2k} terms are given by performing this same process on the sequence  y_n

Method 2: A real DFT and a DCT-4

This method is probably not particularly useful given how the next method works - but I thought it was interesting so I'm listing it. Starting from the beginning again:

 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) k }{N} } \right\}

Rearrange the above into two halves by interleaving the input:

 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) k }{N} } + x_{N-2n-1} \mathrm{e}^{ - \frac{ \mathrm{j} \pi \left( 2n + \frac{1}{2} \right) k }{N }} \left( -1 \right)^k \right\}

Again, we conjugate terms where we know it won't impact the result:

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

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

For the  X_{2k} case, we can write this as:

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

So for the even output indices, we can use a DFT on some rearranged data followed by a twiddle. For the  X_{2k+1} outputs, this particular form gets ugly - but recalling that in the first method, the  X_{2k+1} outputs were not recursive and simply based on a DCT-4, we just steal that result:

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

I won't bother explicitly listing the steps for this one. They are simple to figure out.

Method 3: A single real DFT

This is almost certainly the method you should use to compute the DCT-2. Recall the third equation in the second method:

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

I didn't spot this immediately, but the summation term also happens to be the first stage of a decimation in frequency transform. Recall the DIF FFT stage:

 Y_k = \displaystyle\sum\limits_{n=0}^{N-1} y_n \mathrm{e}^{ - \frac{\mathrm{j} 2 \pi n k }{N} }

 Y_k = \displaystyle\sum\limits_{n=0}^{N/2-1} y_n \mathrm{e}^{ - \frac{\mathrm{j} 2 \pi n k }{N} } + y_{n+N/2} \mathrm{e}^{ - \frac{\mathrm{j} 2 \pi n k }{N} } \left( -1 \right)^k

So we create a new vector  y_n :

 y_n = \left\{ \begin{array}{l l} x_{2n} & \quad {0 <= n < N/2} \\ x_{N-2n-1} & \quad {N/2 <= n < N} \end{array} \right.

Then substitute back into  X_k :

 X_k = \Re \left\{ \mathrm{e}^{ - \frac{ \mathrm{j} \pi k }{ 2 N } } \displaystyle\sum\limits_{n=0}^{N-1} y_n \mathrm{e}^{ - \frac{ \mathrm{j} 2 \pi n k }{N} } \right\}

This actually becomes an  N/2 length DFT if you use a real transform.

Derivation of fast DCT-4 algorithm based on DFT

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.