Discrete cosine transform

Trigonometric relation

The relation

\[\ctwiddle{ 2 \pi }{ \left( m + \frac{1}{2} \right) \left( 2 k + l \right) }{ 2 N } = \left( - 1 \right)^l \ctwiddle{ 2 \pi }{ \left( n + \frac{1}{2} \right) \left( 2 k + l \right) }{ 2 N }\]

where \(m \equiv N - 1 - n\) is repeatedly used in this section, where the symbols \(n, k, l\) are all integers.

Derivation
\[ \begin{align}\begin{aligned}\ctwiddle{ 2 \pi }{ \left( m + \frac{1}{2} \right) \left( 2 k + l \right) }{ 2 N } & = \ctwiddle{ 2 \pi }{ \left\{ \left( N - 1 - n \right) + \frac{1}{2} \right\} \left( 2 k + l \right) }{ 2 N }\\& = \ctwiddle{ 2 \pi }{ \left\{ N - \left( n + \frac{1}{2} \right) \right\} \left( 2 k + l \right) }{ 2 N }\\& = \cos \left( \left( 2 k + l \right) \pi - 2 \pi \frac{ \left( n + \frac{1}{2} \right) \left( 2 k + l \right) }{ 2 N } \right).\end{aligned}\end{align} \]

By utilizing

\[ \begin{align}\begin{aligned}\cos \left( \alpha - \beta \right) & = \cos \alpha \cos \beta + \sin \alpha \sin \beta,\\\cos \left( \left( 2 k + l \right) \pi \right) & = \left( - 1 \right)^l,\\\sin \left( \left( 2 k + l \right) \pi \right) & = 0,\end{aligned}\end{align} \]

we end up with

\[\ctwiddle{ 2 \pi }{ \left( m + \frac{1}{2} \right) \left( 2 k + l \right) }{ 2 N } = \left( - 1 \right)^l \ctwiddle{ 2 \pi }{ \left( n + \frac{1}{2} \right) \left( 2 k + l \right) }{ 2 N }.\]

Forward transform

We consider the forward transform (discrete cosine transform of type II):

\[X_k \equiv 2 \sum_{n = 0}^{N - 1} x_n \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 N} \equiv \dctii{k}{N}{x_0}{x_1}{x_{N - 1}}\]

for \(\seq{k}{0}{1}{N - 1}\).

Assuming that \(N\) is a multiple of \(2\), we consider the even and odd wavenumbers separately:

\[ \begin{align}\begin{aligned}X_{2 k} = & 2 \sum_{n = 0}^{N - 1} x_n \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k \right)}{2 N},\\X_{2 k + 1} = & 2 \sum_{n = 0}^{N - 1} x_n \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N},\end{aligned}\end{align} \]

for \(\seq{k}{0}{1}{N / 2 - 1}\).

The summation with respect to \(n\) is split into two parts as well following

\[ \begin{align}\begin{aligned}\sum_{n = 0}^{N - 1} q_n = & \sum_{n = 0}^{N / 2 - 1} q_n + \sum_{n = N / 2}^{N - 1} q_n,\\= & \sum_{n = 0}^{N / 2 - 1} q_n + \sum_{n = 0}^{N / 2 - 1} q_{N - 1 - n},\end{aligned}\end{align} \]

giving

\[ \begin{align}\begin{aligned}X_{2 k} = & 2 \sum_{n = 0}^{N / 2 - 1} x_n \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k \right)}{2 N} + 2 \sum_{n = 0}^{N / 2 - 1} x_{N - 1 - n} \ctwiddle{2 \pi}{\left( N - 1 - n + \frac{1}{2} \right) \left( 2 k \right)}{2 N}\\= & 2 \sum_{n = 0}^{N / 2 - 1} x_n \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k \right)}{2 N} + 2 \sum_{n = 0}^{N / 2 - 1} x_{N - 1 - n} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k \right)}{2 N}\\= & 2 \sum_{n = 0}^{N / 2 - 1} \left( x_n + x_{N - 1 - n} \right) \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 \left( N / 2 \right)}\\= & \dctii{k}{N / 2}{x_0 + x_{N - 1}}{x_1 + x_{N - 2}}{x_{N / 2 - 1} + x_{N / 2}},\end{aligned}\end{align} \]

and

\[ \begin{align}\begin{aligned}X_{2 k + 1} = & 2 \sum_{n = 0}^{N / 2 - 1} x_n \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N} + 2 \sum_{n = 0}^{N / 2 - 1} x_{N - 1 - n} \ctwiddle{2 \pi}{\left( N - 1 - n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\\= & 2 \sum_{n = 0}^{N / 2 - 1} x_n \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N} - 2 \sum_{n = 0}^{N / 2 - 1} x_{N - 1 - n} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\\= & 2 \sum_{n = 0}^{N / 2 - 1} \left( x_n - x_{N - 1 - n} \right) \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N},\end{aligned}\end{align} \]

for \(\seq{k}{0}{1}{N / 2 - 1}\). Note that we adopt the trigonometric relations derived before.

By utilizing the product-to-sum identity:

\[\cos \alpha \equiv \frac{ \cos \left( \alpha + \beta \right) + \cos \left( \alpha - \beta \right) }{ 2 \cos \beta }\]

with

\[\beta = 2 \pi \frac{ n + \frac{1}{2} }{ 2 N },\]

we obtain

\[ \begin{align}\begin{aligned}X_{2 k + 1} = & 2 \sum_{n = 0}^{N / 2 - 1} \frac{ x_n - x_{N - 1 - n} }{2 \cos \beta} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( k + 1 \right)}{2 \left( N / 2 \right)} + 2 \sum_{n = 0}^{N / 2 - 1} \frac{ x_n - x_{N - 1 - n} }{2 \cos \beta} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 \left( N / 2 \right)}\\= & \dctii{ k + 1 }{ N / 2 }{ \frac{ x_0 - x_{N - 1} }{2 \cos \beta} }{ \frac{ x_1 - x_{N - 2} }{2 \cos \beta} }{ \frac{ x_{N / 2 - 1} - x_{N / 2} }{2 \cos \beta} } + \dctii{ k }{ N / 2 }{ \frac{ x_0 - x_{N - 1} }{2 \cos \beta} }{ \frac{ x_1 - x_{N - 2} }{2 \cos \beta} }{ \frac{ x_{N / 2 - 1} - x_{N / 2} }{2 \cos \beta} },\end{aligned}\end{align} \]

which can be applied directly for \(\seq{k}{0}{1}{N / 2 - 2}\). For \(k = N / 2 - 1\) (i.e. a corner case), assigning \(k = N / 2 - 1\) to the first term reveals

\[ \begin{align}\begin{aligned}\dctii{ N / 2 }{ N / 2 }{ \frac{ x_0 - x_{N - 1} }{2 \cos \beta} }{ \frac{ x_1 - x_{N - 2} }{2 \cos \beta} }{ \frac{ x_{N / 2 - 1} - x_{N / 2} }{2 \cos \beta} } & = \sum_{n = 0}^{N / 2 - 1} \frac{ x_n - x_{N - 1 - n} }{2 \cos \beta} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) N / 2}{2 \left( N / 2 \right)}\\& = \sum_{n = 0}^{N / 2 - 1} \frac{ x_n - x_{N - 1 - n} }{2 \cos \beta} \cos \left( \frac{2 n + 1}{2} \pi \right)\\& = 0.\end{aligned}\end{align} \]

In summary, we obtain the following recurrence relation:

\[ \begin{align}\begin{aligned}X_{2 k} & = \dctii{k}{N / 2}{x_{0} + x_{N - 1}}{x_{1} + x_{N - 2}}{x_{N / 2 - 1} + x_{N / 2}},\\X_{2 k + 1} & = \dctii{ k }{ N / 2 }{ \frac{ x_0 - x_{N - 1} }{2 \cos \beta} }{ \frac{ x_1 - x_{N - 2} }{2 \cos \beta} }{ \frac{ x_{N / 2 - 1} - x_{N / 2} }{2 \cos \beta} }\\& + \dctii{ k + 1 }{ N / 2 }{ \frac{ x_0 - x_{N - 1} }{2 \cos \beta} }{ \frac{ x_1 - x_{N - 2} }{2 \cos \beta} }{ \frac{ x_{N / 2 - 1} - x_{N / 2} }{2 \cos \beta} },\end{aligned}\end{align} \]

which is valid for \(\seq{k}{0}{1}{N / 2 - 1}\). Recall that

\[\dctii{ N / 2 }{ N / 2 }{ \frac{ x_0 - x_{N - 1} }{2 \cos \beta} }{ \frac{ x_1 - x_{N - 2} }{2 \cos \beta} }{ \frac{ x_{N / 2 - 1} - x_{N / 2} }{2 \cos \beta} } = 0.\]

Note that, \(X_0 = 2 x_0\) for \(N = 1\).

src/dct.c
37static int dct2 (
38    const size_t nitems,
39    const size_t stride,
40    const double * const restrict table,
41    double * const restrict xs,
42    double * const restrict ys
43) {
44  if (1 == nitems) {
45    xs[0] *= 2.;
46    return 0;
47  }
48  if (0 == nitems % 2) {
49    // divide and conquer
50    const size_t nh = nitems / 2;
51    double * const buf0 = ys +  0;
52    double * const buf1 = ys + nh;
53    // create input buffers of DCT-II
54    for (size_t n = 0; n < nh; n++) {
55      // c: 1 / [ 2 cos(beta) ]
56      const double c = table[(2 * n + 1) * stride];
57      const double value0 = xs[             n];
58      const double value1 = xs[nitems - 1 - n];
59      buf0[n] = 1. * (value0 + value1);
60      buf1[n] = c  * (value0 - value1);
61    }
62    // solve two sub problems
63    dct2(nh, stride * 2, table, buf0, xs);
64    dct2(nh, stride * 2, table, buf1, xs);
65    // distribute results
66    for (size_t k = 0; k < nh; k++) {
67      // to even frequencies
68      xs[k * 2 + 0] = buf0[k];
69      // to odd frequencies
70      // for the last element, "k+1"-th element is zero
71      const double value0 =                    buf1[k    ];
72      const double value1 = nh - 1 == k ? 0. : buf1[k + 1];
73      xs[k * 2 + 1] = value0 + value1;
74    }
75  } else {
76    // fallback to N^2 DCT2
77    for (size_t k = 0; k < nitems; k++) {
78      double * const y = ys + k;
79      *y = 0.;
80      for (size_t n = 0; n < nitems; n++) {
81        const double phase = pi * (2 * n + 1) * k / (2 * nitems);
82        *y += 2. * xs[n] * cos(phase);
83      }
84    }
85    for (size_t k = 0; k < nitems; k++) {
86      xs[k] = ys[k];
87    }
88  }
89  return 0;
90}

Backward transform

We consider the inverse transform (discrete cosine transform of type III):

\[x_n = \frac{1}{N} \left\{ \frac{\hat{X}_0}{2} + \sum_{k = 1}^{N - 1} \hat{X}_k \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 N} \right\}\]

for \(\seq{n}{0}{1}{N - 1}\). For later convenience, we define

\[\begin{split}X_k = \begin{cases} \frac{\hat{X}_0}{2} & k = 0, \\ \hat{X}_k & \text{otherwise}, \end{cases}\end{split}\]

to write the inverse transform as

\[x_n = \frac{1}{N} \sum_{k = 0}^{N - 1} X_k \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 N} \equiv \dctiii{n}{N}{X_0}{X_1}{X_{N - 1}}\]

for \(\seq{n}{0}{1}{N - 1}\).

Note that we perform this conversion a priori:

src/dct.c
215xs[0] *= 0.5;
216dct3(nitems, 1, table, xs, ys);

Assuming that \(N\) is a multiple of \(2\), we decompose the right-hand side into two components:

\[ \begin{align}\begin{aligned}x_n & = \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k \right)}{2 N} + \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\\& = \frac{1}{2} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} X_{2 k} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 \left( N / 2 \right)} + \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\\& = \frac{1}{2} \dctiii{n}{N / 2}{X_0}{X_2}{X_{N - 2}} + \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}.\end{aligned}\end{align} \]

To decompose \(n\) as well, we consider \(n \leftarrow N - 1 - n\) to yield

\[ \begin{align}\begin{aligned}x_{N - 1 - n} & = \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k} \ctwiddle{2 \pi}{\left( N - n - \frac{1}{2} \right) \left( 2 k \right)}{2 N} + \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( N - n - \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\\& = \frac{1}{2} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} X_{2 k} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 \left( N / 2 \right)} - \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\\& = \frac{1}{2} \dctiii{n}{N / 2}{X_0}{X_2}{X_{N - 2}} - \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N},\end{aligned}\end{align} \]

due to one of the trigonometric relations derived before. By adopting the product-to-sum identity:

\[\cos \alpha \equiv \frac{ \cos \left( \alpha + \beta \right) + \cos \left( \alpha - \beta \right) }{ 2 \cos \beta }\]

with

\[\beta = 2 \pi \frac{ n + \frac{1}{2} }{ 2 N },\]

we obtain

\[ \begin{align}\begin{aligned}& \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\\= & \frac{1}{2 N \cos \beta} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k \right)}{2 N} + \frac{1}{2 N \cos \beta} \sum_{k = 1}^{N / 2} X_{2 k - 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k \right)}{2 N}.\end{aligned}\end{align} \]

The first term leads to

\[\frac{1}{4 \cos \beta} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 \left( N / 2 \right)} = \frac{1}{4 \cos \beta} \dctiii{n}{N / 2}{X_1}{X_3}{X_{N - 1}},\]

while the second term is

\[ \begin{align}\begin{aligned}& \frac{1}{4 \cos \beta} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} X_{2 k - 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 \left( N / 2 \right)} - \frac{1}{2 N \cos \beta} X_{-1} + \frac{1}{2 N \cos \beta} X_{N - 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) N}{2 N}\\= & \frac{1}{4 \cos \beta} \dctiii{n}{N / 2}{X_{-1}}{X_1}{X_{N - 3}},\end{aligned}\end{align} \]

where we artificially specify \(X_{-1} = 0\). As a consequence, we obtain

\[\frac{1}{N} \sum_{k = 0}^{N / 2 - 1} X_{2 k + 1} \ctwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N} = \frac{1}{4 \cos \beta} \dctiii{n}{N / 2}{X_{-1} + X_1}{X_1 + X_3}{X_{N - 3} + X_{N - 1}}.\]

In summary, we obtain the following recurrence relation:

\[ \begin{align}\begin{aligned}x_n & = \frac{1}{2} \dctiii{n}{N / 2}{X_0}{X_2}{X_{N - 2}} + \frac{1}{4 \cos \beta} \dctiii{n}{N / 2}{X_{-1} + X_1}{X_1 + X_3}{X_{N - 3} + X_{N - 1}},\\x_{N - 1 - n} & = \frac{1}{2} \dctiii{n}{N / 2}{X_0}{X_2}{X_{N - 2}} - \frac{1}{4 \cos \beta} \dctiii{n}{N / 2}{X_{-1} + X_1}{X_1 + X_3}{X_{N - 3} + X_{N - 1}},\end{aligned}\end{align} \]

for \(\seq{n}{0}{1}{N / 2 - 1}\). Note that \(x_0 = X_0\) for \(N = 1\).

src/dct.c
 93static int dct3 (
 94    const size_t nitems,
 95    const size_t stride,
 96    const double * const restrict table,
 97    double * const restrict xs,
 98    double * const restrict ys
 99) {
100  if (1 == nitems) {
101    return 0;
102  }
103  if (0 == nitems % 2) {
104    // divide and conquer
105    const size_t nh = nitems / 2;
106    double * const buf0 = ys +  0;
107    double * const buf1 = ys + nh;
108    // create input buffers of DCT-III
109    for (size_t k = 0; k < nh; k++) {
110      buf0[k] = xs[k * 2];
111      const double value0 = 0 == k ? 0. : xs[k * 2 - 1];
112      const double value1 =               xs[k * 2 + 1];
113      buf1[k] = value0 + value1;
114    }
115    // solve two sub problems
116    dct3(nh, stride * 2, table, buf0, xs);
117    dct3(nh, stride * 2, table, buf1, xs);
118    // combine results of sub problems
119    for (size_t n = 0; n < nh; n++) {
120      // c: 1 / [ 2 cos(beta) ]
121      const double c = table[(2 * n + 1) * stride];
122      const double value0 = 0.5 * buf0[n];
123      const double value1 = 0.5 * c * buf1[n];
124      xs[             n] = value0 + value1;
125      xs[nitems - 1 - n] = value0 - value1;
126    }
127  } else {
128    // fallback to N^2 DCT3
129    for (size_t n = 0; n < nitems; n++) {
130      double * const y = ys + n;
131      *y = 0.;
132      for (size_t k = 0; k < nitems; k++) {
133        const double phase = pi * (2 * n + 1) * k / (2 * nitems);
134        *y += xs[k] * cos(phase);
135      }
136      *y /= 1. * nitems;
137    }
138    for (size_t n = 0; n < nitems; n++) {
139      xs[n] = ys[n];
140    }
141  }
142  return 0;
143}

Reference

  • Lee, “A New Algorithm to Compute the Discrete Cosine Transform”, IEEE T. Acoust. Speech, 1984