Discrete cosine transform

The DCT is widely used in many applications, which include the following.

The DCT-II, also known as simply the DCT, is the most important image compression technique. It is used in image compression standards such as JPEG, and video compression standards such as H.26x, MJPEG, MPEG, DV, Theora and Daala. There, the two-dimensional DCT-II of ${\displaystyle N\times N}$ blocks are computed and the results are quantized and entropy coded. In this case, ${\displaystyle N}$ is typically 8 and the DCT-II formula is applied to each row and column of the block. The result is an 8 × 8 transform coefficient array in which the ${\displaystyle (0,0)}$ element (top-left) is the DC (zero-frequency) component and entries with increasing vertical and horizontal index values represent higher vertical and horizontal spatial frequencies.

Advanced Video Coding (AVC) uses the integer DCT[23][1] (IntDCT), an integer approximation of the DCT.[2][1] It uses 4x4 and 8x8 integer DCT blocks. High Efficiency Video Coding (HEVC) and the High Efficiency Image Format (HEIF) use varied integer DCT block sizes between 4x4 and 32x32 pixels.[4][1] As of 2019, AVC is by far the most commonly used format for the recording, compression and distribution of video content, used by 91% of video developers, followed by HEVC which is used by 43% of developers.[54]

Multidimensional DCTs (MD DCTs) have several applications, mainly 3-D DCTs such as the 3-D DCT-II, which has several new applications like Hyperspectral Imaging coding systems,[92] variable temporal length 3-D DCT coding,[93] video coding algorithms,[94] adaptive video coding [95] and 3-D Compression.[96] Due to enhancement in the hardware, software and introduction of several fast algorithms, the necessity of using M-D DCTs is rapidly increasing. DCT-IV has gained popularity for its applications in fast implementation of real-valued polyphase filtering banks,[97] lapped orthogonal transform[98][99] and cosine-modulated wavelet bases.[100]

DCT plays a very important role in digital signal processing. By using the DCT, the signals can be compressed. DCT can be used in electrocardiography for the compression of ECG signals. DCT2 provides a better compression ratio than DCT.

The DCT is widely implemented in digital signal processors (DSP), as well as digital signal processing software. Many companies have developed DSPs based on DCT technology. DCTs are widely used for applications such as encoding, decoding, video, audio, multiplexing, control signals, signaling, and analog-to-digital conversion. DCTs are also commonly used for high-definition television (HDTV) encoder/decoder chips.[1]

A common issue with DCT compression in digital media are blocky compression artifacts,[101] caused by DCT blocks.[3] The DCT algorithm can cause block-based artifacts when heavy compression is applied. Due to the DCT being used in the majority of digital image and video coding standards (such as the JPEG, H.26x and MPEG formats), DCT-based blocky compression artifacts are widespread in digital media. In a DCT algorithm, an image (or frame in an image sequence) is divided into square blocks which are processed independently from each other, then the DCT of these blocks is taken, and the resulting DCT coefficients are quantized. This process can cause blocking artifacts, primarily at high data compression ratios.[101] This can also cause the "mosquito noise" effect, commonly found in digital video (such as the MPEG formats).[102]

DCT blocks are often used in glitch art.[3] The artist Rosa Menkman makes use of DCT-based compression artifacts in her glitch art,[103] particularly the DCT blocks found in most digital media formats such as JPEG digital images and MP3 digital audio.[3] Another example is Jpegs by German photographer Thomas Ruff, which uses intentional JPEG artifacts as the basis of the picture's style.[104][105]

${\displaystyle X_{k}={\frac {1}{2}}(x_{0}+(-1)^{k}x_{N-1})+\sum _{n=1}^{N-2}x_{n}\cos \left[{\frac {\pi }{N-1}}nk\right]\quad \quad k=0,\dots ,N-1.}$

Some authors further multiply the x₀ and x_N-1 terms by √2, and correspondingly multiply the X₀ and X_N-1 terms by 1/√2. This makes the DCT-I matrix orthogonal, if one further multiplies by an overall scale factor of ${\displaystyle {\sqrt {\tfrac {2}{N-1}}}}$, but breaks the direct correspondence with a real-even DFT.

The DCT-I is exactly equivalent (up to an overall scale factor of 2), to a DFT of ${\displaystyle 2N-2}$ real numbers with even symmetry. For example, a DCT-I of N=5 real numbers abcde is exactly equivalent to a DFT of eight real numbers abcdedcb (even symmetry), divided by two. (In contrast, DCT types II-IV involve a half-sample shift in the equivalent DFT.)

Note, however, that the DCT-I is not defined for N less than 2. (All other DCT types are defined for any positive N.)

Thus, the DCT-I corresponds to the boundary conditions: x_n is even around n = 0 and even around n = N−1; similarly for X_k.

${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}\cos \left[{\frac {\pi }{N}}\left(n+{\frac {1}{2}}\right)k\right]\quad \quad k=0,\dots ,N-1.}$

The DCT-II is probably the most commonly used form, and is often simply referred to as "the DCT".[5][6]

This transform is exactly equivalent (up to an overall scale factor of 2) to a DFT of ${\displaystyle 4N}$ real inputs of even symmetry where the even-indexed elements are zero. That is, it is half of the DFT of the ${\displaystyle 4N}$ inputs ${\displaystyle y_{n}}$, where ${\displaystyle y_{2n}=0}$, ${\displaystyle y_{2n+1}=x_{n}}$ for ${\displaystyle 0\leq n<N}$, ${\displaystyle y_{2N}=0}$, and ${\displaystyle y_{4N-n}=y_{n}}$ for ${\displaystyle 0<n<2N}$. DCT II transformation is also possible using 2N signal followed by a multiplication by half shift. This is demonstrated by Makhoul.

Some authors further multiply the X₀ term by 1/√2 and multiply the resulting matrix by an overall scale factor of ${\displaystyle {\sqrt {\tfrac {2}{N}}}}$ (see below for the corresponding change in DCT-III). This makes the DCT-II matrix orthogonal, but breaks the direct correspondence with a real-even DFT of half-shifted input. This is the normalization used by Matlab, for example.[106] In many applications, such as JPEG, the scaling is arbitrary because scale factors can be combined with a subsequent computational step (e.g. the quantization step in JPEG[107]), and a scaling can be chosen that allows the DCT to be computed with fewer multiplications.[108][109]

The DCT-II implies the boundary conditions: x_n is even around n = −1/2 and even around n = N−1/2; X_k is even around k = 0 and odd around k = N.

${\displaystyle X_{k}={\frac {1}{2}}x_{0}+\sum _{n=1}^{N-1}x_{n}\cos \left[{\frac {\pi }{N}}n\left(k+{\frac {1}{2}}\right)\right]\quad \quad k=0,\dots ,N-1.}$

Because it is the inverse of DCT-II (up to a scale factor, see below), this form is sometimes simply referred to as "the inverse DCT" ("IDCT").[6]

Some authors divide the x₀ term by √2 instead of by 2 (resulting in an overall x₀/√2 term) and multiply the resulting matrix by an overall scale factor of ${\displaystyle {\sqrt {\tfrac {2}{N}}}}$ (see above for the corresponding change in DCT-II), so that the DCT-II and DCT-III are transposes of one another. This makes the DCT-III matrix orthogonal, but breaks the direct correspondence with a real-even DFT of half-shifted output.

The DCT-III implies the boundary conditions: x_n is even around n = 0 and odd around n = N; X_k is even around k = −1/2 and even around k = N−1/2.

${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}\cos \left[{\frac {\pi }{N}}\left(n+{\frac {1}{2}}\right)\left(k+{\frac {1}{2}}\right)\right]\quad \quad k=0,\dots ,N-1.}$

The DCT-IV matrix becomes orthogonal (and thus, being clearly symmetric, its own inverse) if one further multiplies by an overall scale factor of ${\displaystyle {\sqrt {2/N}}}$.

A variant of the DCT-IV, where data from different transforms are overlapped, is called the modified discrete cosine transform (MDCT).[110]

The DCT-IV implies the boundary conditions: x_n is even around n = −1/2 and odd around n = N−1/2; similarly for X_k.

DCTs of types I-IV treat both boundaries consistently regarding the point of symmetry: they are even/odd around either a data point for both boundaries or halfway between two data points for both boundaries. By contrast, DCTs of types V-VIII imply boundaries that are even/odd around a data point for one boundary and halfway between two data points for the other boundary.

In other words, DCT types I-IV are equivalent to real-even DFTs of even order (regardless of whether N is even or odd), since the corresponding DFT is of length 2(N−1) (for DCT-I) or 4N (for DCT-II/III) or 8N (for DCT-IV). The four additional types of discrete cosine transform[111] correspond essentially to real-even DFTs of logically odd order, which have factors of N ± ½ in the denominators of the cosine arguments.

However, these variants seem to be rarely used in practice. One reason, perhaps, is that FFT algorithms for odd-length DFTs are generally more complicated than FFT algorithms for even-length DFTs (e.g. the simplest radix-2 algorithms are only for even lengths), and this increased intricacy carries over to the DCTs as described below.

(The trivial real-even array, a length-one DFT (odd length) of a single number a, corresponds to a DCT-V of length N = 1.)

For example, a two-dimensional DCT-II of an image or a matrix is simply the one-dimensional DCT-II, from above, performed along the rows and then along the columns (or vice versa). That is, the 2D DCT-II is given by the formula (omitting normalization and other scale factors, as above):

${\displaystyle {\begin{aligned}X_{k_{1},k_{2}}&=\sum _{n_{1}=0}^{N_{1}-1}\left(\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1},n_{2}}\cos \left[{\frac {\pi }{N_{2}}}\left(n_{2}+{\frac {1}{2}}\right)k_{2}\right]\right)\cos \left[{\frac {\pi }{N_{1}}}\left(n_{1}+{\frac {1}{2}}\right)k_{1}\right]\\&=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1},n_{2}}\cos \left[{\frac {\pi }{N_{1}}}\left(n_{1}+{\frac {1}{2}}\right)k_{1}\right]\cos \left[{\frac {\pi }{N_{2}}}\left(n_{2}+{\frac {1}{2}}\right)k_{2}\right].\end{aligned}}}$

The inverse of a multi-dimensional DCT is just a separable product of the inverses of the corresponding one-dimensional DCTs (see above), e.g. the one-dimensional inverses applied along one dimension at a time in a row-column algorithm.

The 3-D DCT-II is only the extension of 2-D DCT-II in three dimensional space and mathematically can be calculated by the formula

${\displaystyle {\begin{aligned}X_{k_{1},k_{2},k_{3}}&=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}\sum _{n_{3}=0}^{N_{3}-1}x_{n_{1},n_{2},n_{3}}\cos \left[{\frac {\pi }{N_{1}}}\left(n_{1}+{\frac {1}{2}}\right)k_{1}\right]\cos \left[{\frac {\pi }{N_{2}}}\left(n_{2}+{\frac {1}{2}}\right)k_{2}\right]\cos \left[{\frac {\pi }{N_{3}}}\left(n_{3}+{\frac {1}{2}}\right)k_{3}\right],\quad \forall k_{i}=0,1,2,\dots ,N_{i}-1.\end{aligned}}}$

The inverse of 3-D DCT-II is 3-D DCT-III and can be computed from the formula given by

${\displaystyle {\begin{aligned}x_{n_{1},n_{2},n_{3}}&=\sum _{k_{1}=0}^{N_{1}-1}\sum _{k_{2}=0}^{N_{2}-1}\sum _{k_{3}=0}^{N_{3}-1}X_{k_{1},k_{2},k_{3}}\cos \left[{\frac {\pi }{N_{1}}}\left(n_{1}+{\frac {1}{2}}\right)k_{1}\right]\cos \left[{\frac {\pi }{N_{2}}}\left(n_{2}+{\frac {1}{2}}\right)k_{2}\right]\cos \left[{\frac {\pi }{N_{3}}}\left(n_{3}+{\frac {1}{2}}\right)k_{3}\right],\quad \forall n_{i}=0,1,2,\dots ,N_{i}-1.\end{aligned}}}$

Technically, computing a two-, three- (or -multi) dimensional DCT by sequences of one-dimensional DCTs along each dimension is known as a row-column algorithm. As with multidimensional FFT algorithms, however, there exist other methods to compute the same thing while performing the computations in a different order (i.e. interleaving/combining the algorithms for the different dimensions). Owing to the rapid growth in the applications based on the 3-D DCT, several fast algorithms are developed for the computation of 3-D DCT-II. Vector-Radix algorithms are applied for computing M-D DCT to reduce the computational complexity and to increase the computational speed. To compute 3-D DCT-II efficiently, a fast algorithm, Vector-Radix Decimation in Frequency (VR DIF) algorithm was developed.

In order to apply the VR DIF algorithm the input data is to be formulated and rearranged as follows.[112][113] The transform size N x N x N is assumed to be 2.

The four basic stages of computing 3-D DCT-II using VR DIF Algorithm.

${\displaystyle {\begin{array}{lcl}{\tilde {x}}(n_{1},n_{2},n_{3})=x(2n_{1},2n_{2},2n_{3})\\{\tilde {x}}(n_{1},n_{2},N-n_{3}-1)=x(2n_{1},2n_{2},2n_{3}+1)\\{\tilde {x}}(n_{1},N-n_{2}-1,n_{3})=x(2n_{1},2n_{2}+1,2n_{3})\\{\tilde {x}}(n_{1},N-n_{2}-1,N-n_{3}-1)=x(2n_{1},2n_{2}+1,2n_{3}+1)\\{\tilde {x}}(N-n_{1}-1,n_{2},n_{3})=x(2n_{1}+1,2n_{2},2n_{3})\\{\tilde {x}}(N-n_{1}-1,n_{2},N-n_{3}-1)=x(2n_{1}+1,2n_{2},2n_{3}+1)\\{\tilde {x}}(N-n_{1}-1,N-n_{2}-1,n_{3})=x(2n_{1}+1,2n_{2}+1,2n_{3})\\{\tilde {x}}(N-n_{1}-1,N-n_{2}-1,N-n_{3}-1)=x(2n_{1}+1,2n_{2}+1,2n_{3}+1)\\\end{array}}}$where ${\displaystyle 0\leq n_{1},n_{2},n_{3}\leq {\frac {N}{2}}-1}$

The figure to the adjacent shows the four stages that are involved in calculating 3-D DCT-II using VR DIF algorithm. The first stage is the 3-D reordering using the index mapping illustrated by the above equations. The second stage is the butterfly calculation. Each butterfly calculates eight points together as shown in the figure just below, where ${\displaystyle c(\phi _{i})=\cos(\phi _{i})}$.

The original 3-D DCT-II now can be written as

${\displaystyle X(k_{1},k_{2},k_{3})=\sum _{n_{1}=1}^{N-1}\sum _{n_{2}=1}^{N-1}\sum _{n_{3}=1}^{N-1}{\tilde {x}}(n_{1},n_{2},n_{3})\cos(\phi k_{1})\cos(\phi k_{2})\cos(\phi k_{3})}$

where ${\displaystyle \phi _{i}={\frac {\pi }{2N}}(4N_{i}+1),{\text{ and }}i=1,2,3}$ .

If the even and the odd parts of ${\displaystyle k_{1},k_{2}}$ and ${\displaystyle k_{3}}$ and are considered, the general formula for the calculation of the 3-D DCT-II can be expressed as

The single butterfly stage of VR DIF algorithm.

${\displaystyle X(k_{1},k_{2},k_{3})=\sum _{n_{1}=1}^{{\tfrac {N}{2}}-1}\sum _{n_{2}=1}^{{\tfrac {N}{2}}-1}\sum _{n_{1}=1}^{{\tfrac {N}{2}}-1}{\tilde {x}}_{ijl}(n_{1},n_{2},n_{3})\cos(\phi (2k_{1}+i)\cos(\phi (2k_{2}+j)\cos(\phi (2k_{3}+l))}$

where

${\displaystyle {\tilde {x}}_{ijl}(n_{1},n_{2},n_{3})={\tilde {x}}(n_{1},n_{2},n_{3})+(-1)^{l}{\tilde {x}}\left(n_{1},n_{2},n_{3}+{\frac {n}{2}}\right)}$

${\displaystyle +(-1)^{j}{\tilde {x}}\left(n_{1},n_{2}+{\frac {n}{2}},n_{3}\right)+(-1)^{j+l}{\tilde {x}}\left(n_{1},n_{2}+{\frac {n}{2}},n_{3}+{\frac {n}{2}}\right)}$

${\displaystyle +(-1)^{i}{\tilde {x}}\left(n_{1}+{\frac {n}{2}},n_{2},n_{3}\right)+(-1)^{i+j}{\tilde {x}}\left(n_{1}+{\frac {n}{2}}+{\frac {n}{2}},n_{2},n_{3}\right)}$

${\displaystyle +(-1)^{i+l}{\tilde {x}}\left(n_{1}+{\frac {n}{2}},n_{2},n_{3}+{\frac {n}{3}}\right)}$

${\displaystyle +(-1)^{i+j+l}{\tilde {x}}\left(n_{1}+{\frac {n}{2}},n_{2}+{\frac {n}{2}},n_{3}+{\frac {n}{2}}\right){\text{ where }}i,j,l=0{\text{ or }}1.}$

The whole 3-D DCT calculation needs ${\displaystyle [\log _{2}N]}$ stages, and each stage involves ${\displaystyle N^{3}/8}$ butterflies. The whole 3-D DCT requires ${\displaystyle \left[(N^{3}/8)\log _{2}N\right]}$ butterflies to be computed. Each butterfly requires seven real multiplications (including trivial multiplications) and 24 real additions (including trivial additions). Therefore, the total number of real multiplications needed for this stage is ${\displaystyle \left[(7/8)N^{3}\log _{2}N\right]}$, and the total number of real additions i.e. including the post-additions (recursive additions) which can be calculated directly after the butterfly stage or after the bit-reverse stage are given by[113] ${\displaystyle \underbrace {\left[{\frac {3}{2}}N^{3}\log _{2}N\right]} _{\text{Real}}+\underbrace {\left[{\frac {3}{2}}N^{3}\log _{2}N-3N^{3}+3N^{2}\right]} _{\text{Recursive}}=\left[{\frac {9}{2}}N^{3}\log _{2}N-3N^{3}+3N^{2}\right]}$.

The conventional method to calculate MD-DCT-II is using a Row-Column-Frame (RCF) approach which is computationally complex and less productive on most advanced recent hardware platforms. The number of multiplications required to compute VR DIF Algorithm when compared to RCF algorithm are quite a few in number. The number of Multiplications and additions involved in RCF approach are given by ${\displaystyle \left[{\frac {3}{2}}N^{3}\log _{2}N\right]}$ and ${\displaystyle \left[{\frac {9}{2}}N^{3}\log _{2}N-3N^{3}+3N^{2}\right]}$ respectively. From Table 1, it can be seen that the total number

of multiplications associated with the 3-D DCT VR algorithm is less than that associated with the RCF approach by more than 40%. In addition, the RCF approach involves matrix transpose and more indexing and data swapping than the new VR algorithm. This makes the 3-D DCT VR algorithm more efficient and better suited for 3-D applications that involve the 3-D DCT-II such as video compression and other 3-D image processing applications. The main consideration in choosing a fast algorithm is to avoid computational and structural complexities. As the technology of computers and DSPs advances, the execution time of arithmetic operations (multiplications and additions) is becoming very fast, and regular computational structure becomes the most important factor.[114] Therefore, although the above proposed 3-D VR algorithm does not achieve the theoretical lower bound on the number of multiplications,[115] it has a simpler computational structure as compared to other 3-D DCT algorithms. It can be implemented in place using a single butterfly and possesses the properties of the Cooley–Tukey FFT algorithm in 3-D. Hence, the 3-D VR presents a good choice for reducing arithmetic operations in the calculation of the 3-D DCT-II while keeping the simple structure that characterize butterfly style Cooley–Tukey FFT algorithms.

Two-dimensional DCT frequencies from the JPEG DCT

The image to the right shows a combination of horizontal and vertical frequencies for an 8 x 8 (${\displaystyle N_{1}=N_{2}=8}$) two-dimensional DCT. Each step from left to right and top to bottom is an increase in frequency by 1/2 cycle. For example, moving right one from the top-left square yields a half-cycle increase in the horizontal frequency. Another move to the right yields two half-cycles. A move down yields two half-cycles horizontally and a half-cycle vertically. The source data (8x8) is transformed to a linear combination of these 64 frequency squares.

The M-D DCT-IV is just an extension of 1-D DCT-IV on to M dimensional domain. The 2-D DCT-IV of a matrix or an image is given by

${\displaystyle X_{k,l}=\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}x_{n,m}\cos \left({\frac {(2m+1)(2k+1)\pi }{4N}}\right)\cos \left({\frac {(2n+1)(2l+1)\pi }{4M}}\right){\text{. where }}k=0,1,2...,N-1{\text{ and }}l=0,1,2...,M-1}$.

We can compute the MD DCT-IV using the regular row-column method or we can use the polynomial transform method[116] for the fast and efficient computation. The main idea of this algorithm is to use the Polynomial Transform to convert the multidimensional DCT into a series of 1-D DCTs directly. MD DCT-IV also has several applications in various fields.