Fast DCT algorithm

Forward DCT

The DCT of a sequence $\{x[m],\;\;\;(m=0,\cdots,N-1)\;\}$ can be implemented by FFT. First we define a new sequence $\{\;y[m],\;\;\;(m=0,\cdots,N-1)\;\}$:

$$\left\{ \begin{array}{ll} y[m]\stackrel{\triangle}{=}x[2m] & ... ...triangle}{=}x[2m+1] & (i=0,\cdots, N/2-1) \end{array} \right.$$

Then the DCT of $x[n]$ can be written as the following (the coefficient $a[n]$ is dropped for now for simplicity):

$$X[n] = \sum_{m=0}^{N-1} x[m] \cos\left( \frac{ (2m+1)n\pi}{2N}\right)$$

$$= \sum_{m=0}^{N/2-1}x[2m] \cos\left( \frac{ (4m+1)n\pi}{2N} \right) + \sum_{m=0}^{N/2-1}x[2m+1] \cos\left( \frac{ (4m+3)n\pi}{2N} \right)$$

$$= \sum_{m=0}^{N/2-1}y[m] \cos\left( \frac{ (4m+1)n\pi}{2N} \right)+ \sum_{m=0}^{N/2-1}y[N-1-m] \cos\left( \frac{ (4m+3)n\pi}{2N}\right)$$

where the first summation is for all even terms and second all odd terms. We define for the second summation $m'\stackrel{\triangle}{=}N-1-m$, then the limits of the summation $0$ and $N/2-1$ for $m$ becomes $N-1$ and $N/2$ for $m'$, and the second summation can be written as

$$\sum_{m'=N/2}^{N-1} y[m'] \cos\left( 2n\pi - \frac{ (4m'+1)... ...m'=N/2}^{N-1} y[m'] \cos\left( \frac{ (4m'+1)n\pi}{2N} \right)$$

Now the two summations in the expression of $X[n]$ can be combined

$$X[n]=\sum_{m=0}^{N-1}y[m] \cos\left( \frac{ (4m+1)n\pi}{2N}\right)$$

Next, consider the DFT of $y[m]$:

$$Y[n]=\sum_{m=0}^{N-1} y[m] e^{-j2\pi mn/N} =\sum_{m=0}^{N-1... ...i mn}{N}\right)-j\;\sin\left(\frac{2\pi mn}{N}\right) \right]$$

If we multiply both sides by

$$e^{-jn\pi/2N}=\cos\left(\frac{n\pi}{2N}\right)-j\;\sin\left(\frac{n\pi}{2N}\right)$$

and take the real part of the result (and keep in mind that both $x[m]$ and $y[m]$ are real), we get:

$$Re[e^{-jn\pi/2N} Y[n]]=\sum_{m=0}^{N-1} y[m] \left[\cos\left... ...sum_{m=0}^{N-1} y[m] \cos\left( \frac{ (4m+1)n\pi}{2N} \right)$$

The last equal sign is due to the trigonometric identity:

$$cos(\alpha+\beta)=cos \alpha \;cos \beta - sin \alpha \;sin \beta$$

This expression for $Re[e^{-jn\pi/2N} Y[n]]$ is identical to that for $X[n]$ above, therefore we get

$$X[n]=Re[\; e^{-jn\pi/2N} Y[n]\;]$$

where $Y[n]$ is the DFT of $y[m]$ (defined from $x[m]$) which can be computed using FFT algorithm with time complexity $O(N log_2 N)$.

In summary, fast forward DCT can be implemented in 3 steps:

• Step 1: Generate a sequence $y[m]$ from the given sequence $x[m]$:
  
  $$\left\{ \begin{array}{ll} y[m]=x[2m] & \\ y[N-1-m]=x[2m+1] & (i=0,\cdots, N/2-1) \end{array} \right.$$
  
  

• step 2: Obtain DFT $Y[n]$ of $y[m]$ using FFT. (As $y[m]$ is real, $Y[n]$ is symmetric and only half of the data points need be computed.)
  
  $$Y[n]={\cal F} [y[m]]$$
  
  

• step 3: Obtain DCT $X[n]$ from $Y[n]$ by
  
  $$X[n]=Re[\; e^{-jn\pi/2N} Y[n]\;]$$
  
  

Inverse DCT

The most obvious way to do inverse DCT is to reverse the order and the mathematical operations of the three steps for the forward DCT:

• step 1: Obtain $Y[n]$ from $X[n]$. In step 3 above there are N equations but 2N variables (both real and imaginary parts of $Y[n]$). However, note that as $y[m]'s$ are real, the real part of its spectrum $Y[n]$ is even (N+1 independent variables) and imaginary part odd (N-1 independent variables). So there are only N variables which can be obtained by solving the N equations.
  
  $$X[n]=Re[\; e^{-jn\pi/2N} Y[n]\;]\;\;\;\;(n=0,\cdots,N-1)$$
  
  

• step 2: Obtain $y[m]$ from $Y[n]$ by inverse DFT also using FFT in $N log_2 N$ complexity.
  
  $$y[m]={\cal F}^{-1} [Y[n]]$$
  
  

• Step 3: Obtain $x[m]$ from $y[m]$ by
  
  $$\left\{ \begin{array}{ll} x[2m]=y[m] & \\ x[2m+1]=y[N-1-m] & (i=0,\cdots, N/2-1) \end{array} \right.$$
  
  

However, there is a more efficient way to do the inverse DCT. Consider first the real part of the inverse DFT of the sequence $X[n]e^{jn\pi/2N}$:

$$Re\;\left[ \sum_{n=0}^{N-1} [X[n]e^{jn\pi/2N}]\;e^{j2\pi mn... ...4m+1)n\pi}{2N}\right) = x[2m]\;\;\;\;\;\;(m=0, \cdots, N-1)$$

This equation gives the inverse DCT of all $N/2$ even data points $x(2m)\;\;\;\;(m=0,\cdots,N/2-1)$. To obtain the odd data points, recall that $x[m]=x[2N-m-1]$, and all odd data points

$$x[2m+1]=x[2N-(2m+1)-1]=x[2(N-m-1)]\;\;\;\;(m=0,\cdots,N/2-1)$$

can be obtained from the second half of the previous equation in reverse order $(m=N-1, N-2, \cdots, N/2)$.

In summary, we have these steps to compute IDCT:

• step 1: Generate a sequence $Y[n]$ from the given DCT sequence $X[n]$:
  
  $$Y[n]=X[n]e^{jn\pi/2N} \;\;\;\;\;(n=0,\cdots,N-1)$$
  
  

• step 2: Obtain $y[m]$ from $Y[n]$ by inverse DFT also using FFT. (Only the real part need be computed.)
  
  $$y[m]=Re[ {\cal F}^{-1} [Y[n]] ]$$
  
  

• Step 3: Obtain $x[m]'s$ from $y[m]'s$ by
  
  $$\left\{ \begin{array}{ll} x[2m]=y[m] & \\ x[2m+1]=y[N-1-m] & (i=0,\cdots, N/2-1) \end{array} \right.$$
  
  

These three steps are mathematically equivalent to the steps of the first method.

See the demos.