Understanding the Modified Discrete Cosine Transform (MDCT)

After playing around with discrete cosine transforms, I thought I would implement an MDCT and document my understanding of how everything works. I use some similar techniques to those used on the Wikipedia page as they are helpful for understanding but will add some scanned drawings which I think help (I'm not even close to being clever enough to get a computer to draw these for me).

Prerequisites

The only real background knowledge which I think is relevant to understanding the MDCT is the data extensions which the DCT-4 transform assumes.

First DCT-4 Basis Function with Shifted 2N Sample Input
First DCT-4 Basis Function with Shifted 2N Sample Input

I'll refer to the above image in the Forward Transform overview, but for the mean time, only pay attention to the solid quarter wave. This is the first basis function (i.e.  k=0 ) for an  N length DCT-4. If the basis is continued past  N , it has a repeating symmetrical pattern (the dashed line in the image) which repeats every  4N . The symmetry is even around  -0.5 and odd around  N-0.5 and holds for every basis function of the DCT-4. i.e. The DCT-4 assumes that the input data continues on forever, repeating itself in the following manner:  x_n ,  -x_{N-n-1} ,  -x_n ,  x_{N-n-1} .

Forward Transform

The MDCT takes in  2N real data points and produces  N real outputs. These inputs are designed to overlap, so the first half of the input data should be the second half of the input data of the previous call. The definition is:

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

It should be trivial to see from the above that the MDCT can be computed using a DCT-4 with an extended number of input data points, all of which have been shifted by half a basis. Go back to the crappy drawing and notice the concatenated  N/2 length sequences  a ,  b ,  c and  d . The total length of this sequence is  2N and begins at  N/2 (or half the length of a basis function). We need to get  b ,  c and  d back into the  N point region if we want to compute the MDCT using a DCT-4, this can be achieved with the following concatenated sequence (I will subscript these sequences with  r to denote a reversal of the sequence):

 - c_r - d , a - b_r

If we take the DCT-4 of this concatenated sequence, we have found the MDCT of the input sequence.

Inverse Transform

The inverse MDCT or IMDCT takes in  N real data points and produces  2N real outputs. In this transform, the outputs should overlap such that the first half of the output should be added to the second half of the output data in the previous call. The definition is:

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

Because we know how the DCT-4 assumes the input and output data repeats in a symmetric pattern, we can get this data trivially in exactly the same fashion as we did in the forward transform. In the following Illustration, we take the output from the forward transform and extend it along the basis:

Extended Projection of the MDCT Output on the First DCT-4 Basis
Extended Projection of the MDCT Output on the First DCT-4 Basis

In output row zero, we can see how to extend the input sequence to obtain the 2N points required. We then see in rows two and three how summing the overlapping blocks causes the aliased sequences to cancel in subsequent calls to the IMDCT.

World's Dumbest C MDCT Implementation

I validated all this actually works with a small C program. Follows are the MDCT/IMDCT implementations I came up with... ignore the "twid" input, I cache the modulation factors for the FFT which gets called in the dct4 routine:

/* state should contain double the number of elements as the input buffer (N)
 * and should have all elements initialized to zero prior to calling. The
 * output buffer is actually the first N elements of state after calling. */
void mdct(double *state, const double *input, double *twid, unsigned lenbits)
{
    unsigned rl = 1u << lenbits;
    unsigned i;
    /* Alias the input data with the previous block. */
    for (i = 0; i < rl / 2; i++) {
        state[i]        = - input[rl/2+i]      - input[rl/2-i-1];
        state[rl/2+i]   =   state[rl+i]        - state[rl+rl-i-1];
    }
    /* Save the input block */
    for (i = 0; i < rl; i++)
        state[rl+i]     = input[i];
    /* DCT-4 */
    dct4(state, lenbits, twid);
}
 
/* state should contain double the number of elements as the input buffer (N)
 * and should have all elements initialized to zero prior to calling. The
 * output buffer is actually the first N elements of state after calling. */
void imdct(double *state, const double *input, double *twid, unsigned lenbits)
{
    unsigned rl = 1u << lenbits;
    unsigned i;
    /* Collect contributions from the previous frame to the output buffer */
    for (i = 0; i < rl / 2; i++) {
        state[i]        = - state[rl+rl/2-i-1];
        state[rl/2+i]   = - state[rl+i];
    }
    /* Load the input and run the DCT-4 */
    for (i = 0; i < rl; i++)
        state[rl+i]     = input[i];
    dct4(state + rl, lenbits, twid);
    /* Sum contributions from this frame to the output buffer and perform the
     * required scaling. */
    for (i = 0; i < rl / 2; i++) {
        state[i]        = (state[i]      + state[rl+rl/2+i]) / rl;
        state[rl/2+i]   = (state[rl/2+i] - state[rl+rl-i-1]) / rl;
    }
}

Windowed MDCT Implementation

Typical MDCT implementations will window the input and output data (this can also be thought of as windowing the basis functions - which I think is a more helpful way to understand what is happening). It is really important to note that the window function must be carefully chosen to ensure that the basis functions remain orthogonal! The window makes the basis functions always begin and end near zero. The process has the side effect of de-normalising the basis functions (unless the window is rectangular) and means there will be a window-dependent scaling factor which will need to be applied at the output to achieve perfect reconstruction. The following images show the second basis function of the MDCT both un-windowed and windowed with a half-sine window (given at the end of the post).

Second MDCT Basis Function
Second MDCT Basis Function
Sine Windowed Second MDCT Basis Function
Sine Windowed Second MDCT Basis Function

In a lossy codec this windowing process is somewhat necessary because if the start and end points are not close to zero, the output is likely to periodically glitch for even the slightest errors in the reconstructed MDCT data. This glitching will occur at the boundaries of the transform (i.e. every N points).

We can work out the necessary conditions for the window to obtain perfect reconstruction using the previous drawings (I'd steer away from equations for this one - it's easier to validate the results visually) by applying a window function split into 4 segments to each of the input blocks. I'll do the generic case for a symmetrical window which is applied to both the input and the output. We split the window (which has a length  2N ) into four segments which will be applied to our original input segments  a ,  b ,  c and  d . Because we are defining this window to be symmetric, we can call the pieces:

 u, v, v_r, u_r

Symmetrical Window Impact on MDCT
Symmetrical Window Impact on MDCT

The above illustration shows how our window segments are applied to the input data and the impact that has on the DCT-4 analysed data blob. Following that is the output segments from two sequential IMDCT calls with the windows applied to the output here as well.

We need to make the overlapping terms equal the required output segment i.e.

 c = v_r \left( d_r u + c v_r \right) + u \left( c u - d_r v_r \right)


 d = u_r \left( d u_r + c_r v \right) + v \left( d v - c_r u_r \right)

It is clear from the above that the necessary condition to achieve reconstruction is for  v_r^2 + u^2 = 1 (which implies in this case that  v^2 + u_r^2 = 1 must also be true).

A simple solution to this is:

 w_n = \sin \frac{ \pi \left( n + 0.5 \right) }{2N}

The output requires a scaling factor of 2 for this window.

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.