Discrete Cosine Transform
(and the JPEG Standard)

M. E. Mauel
Columbia University

Introduction

The discrete Cosine transfrom (DCT) is commonly used in digital signal processing. It converts a sequence of real numbers to a sequence of real numbers of the same length. (Recall, the DFT converts a sequence of real numbers (or complex) to a sequence of complex numbers.) The use of the DCT for image compression was introduced in 1974, and it plays an important role in JPEG and MPEG compression schemes.

In this notebook, the DCT and the processing of two dimensional graphics is examined. In particular, we explore

I've used two references. First is "Image Compression Using the Discrete Cosine Transform" by Andrew Watson, Mathematica Journal 4 (1994) p. 81. This is a wonderful collection of Mathematica routines from which I have based this lecture. Second is "The JPEG Still Picture Compression Standard" by Gregory Wallace, Communications of the ACM (1991) p. 35. Both articles are available online on our class website; however, you do not need to read these articles in order to understand this notebook.

DSP Hardware

Today, the DCT and other digital signal processing algorithms are encoded into special hardware that provides high-speed signal processing at low cost.

A parallel video signal processing chip from Oxford Micro Devices (2002) with headquarters in Monroe, CT.  (See http://www.oxfordmicrodevices.com/.) It is designed for fast JPEG, JPEG-2000, MPEG-2 (DVD) and MPEG-4 (streaming video) compression and decompression.

DCT Matrix

As we learned for the Fourier transform, a good way to begin to understand the DCT is to examine the orthogonal modes used to expand digital data.

dctMatrix[nSize_:8] := Table[If[m == 0, Sqrt[1/2],1] Sqrt[2/nSize] Cos[(2 n +1) m π/(2 nSize)], {m, 0, nSize-1}, {n, 0, nSize-1}]

(dct4 = dctMatrix[4]) // TableForm

This transform is orthogonal as shown below…

dct4  .  Transpose[dct4] // Simplify // TableForm

The DCT of a signal s, of length 4, is

dctMatrix[4] . {s[1], s[2], s[3], s[4]} // TableForm

A simple-to-write DCT…

dct1[x_List] := Module[{mDCT,n},
n = Length[x];
mDCT = N[dctMatrix[n]];
mDCT . x];

dct1[{s[1], s[2], s[3], s[4]}]

dct1[{1,2,3,4,4,3,2,1}]

To create an inverse transform, we simply Transpose the conversion matrix. (The transpose is the inverse since the matrix is orthogonal. Since the DCT is real it is also unitary as is the FFT,)

idct1[x_List] := Module[{mDCT,n},
n = Length[x];
mDCT = Transpose[N[dctMatrix[n]]];
mDCT . x];

idct1[dct1[{s[1], s[2], s[3], s[4]}]] // Simplify // Chop

Modes of Length 8

Show[GraphicsGrid[Partition[(BarChart[#1,DisplayFunction->Identity,PlotRange->{-0.5`,0.5`}]&)/@dctMatrix[8],2]],PlotLabel->"DCT Modes of Length 8"]

dct1[{1,2,3,4,4,3,2,1}]

Example 1

sData=Table[((Sin[((i π 6)/32)])/(1+Abs[i-16])),{i,0,32-1}];
ListPlot[sData,PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"Test Signal Length=32"]

sdct=dct1[N[sData]];
ListPlot[sdct,PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"DCT of Test Signal"]

st=idct1[sdct];
ListPlot[st,PlotRange->All,PlotStyle->PointSize[0.02`]]

sf=Fourier[sData];
a1=ListPlot[((Arg[sf])/°),PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"Phase",DisplayFunction->Identity];
a2=ListPlot[Abs[sf],PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"Magnitude",DisplayFunction->Identity];
a3=ListPlot[Re[sf],PlotRange->All,PlotStyle->{RGBColor[0,0,1],PointSize[0.02`]},PlotLabel->"Real",DisplayFunction->Identity];
a4=ListPlot[Im[sf],PlotRange->All,PlotStyle->{RGBColor[1,0,0],PointSize[0.02`]},PlotLabel->"Imaginary",DisplayFunction->Identity];
Show[GraphicsGrid[{{a1,a2},{a3,a4}}],PlotLabel->"FFT (Complex) of Test Signal"]

DCT via the FFT

Since the Cosine is related to the complex exponetial modes of the discrete Fourier transform, the DCT can be computed using the FFT.

Clear[s];
s1 = Table[s[i], {i, 1, 8}]

{xeven, xodd} = Transpose[Partition[s1, 2]]

xeven ~ Join ~ Reverse[xodd]

xfirst = Take[s1,4];
xlast = Take[s1,-4];
Flatten[Transpose[{xfirst, Reverse[xlast]}]]

dctFourier[x_List] := Module[{n,e0, tFactors, xeven, xodd},
n = Length[x];
{xeven, xodd} = Transpose[Partition[x, 2]];
e0 = Exp[-I π/(2 n)];
tFactors = {1} ~ Join~ Table[Sqrt[2] e0^i, {i, 1, n-1}];
Re[tFactors InverseFourier[Join[xeven, Reverse[xodd]]]]]

idctFourier[x_List] :=  Module[{n,e0,tFactors, xt, xfirst, xlast},
n = Length[x];
e0 = Exp[-I π/(2 n)];
tFactors = {1} ~ Join~ Table[Sqrt[2] e0^i, {i, 1, n-1}];
tFactors = Conjugate[N[tFactors]];
xt = Re[Fourier[tFactors x]];
xfirst = Take[xt, n/2];
xlast = Take[xt,-n/2];
Flatten[Transpose[{xfirst, Reverse[xlast]}]]]

sdct2=dctFourier[sData];
ListPlot[sdct2,PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"DCT of Test Signal"]

st2=idctFourier[sdct2];
ListPlot[st2,PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"iDCT of Test Signal"]

sdct - sdct2 // Chop

Sample

nSize=(2^10);
s2=Table[(Sin[(((i-1) 2 π 6)/nSize)]+Sin[(((i-1) 2 π 15)/nSize)]) Sin[((π i)/nSize)],{i,1,nSize}];
p1=ListPlot[s2,PlotLabel->"1k Test Signal"]

sf=Fourier[s2];
a1=ListPlot[((Arg[sf])/°),PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"Phase",DisplayFunction->Identity];
a2=ListPlot[Abs[sf],PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"Magnitude",DisplayFunction->Identity];
Show[GraphicsRow[{a1,a2}]]

ListPlot[Abs[sf[[Range[1,50]]]],PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"FFT"]

ListPlotComplex[sf[[Range[1,50]]], PlotRange-> All, PlotStyle-> PointSize[0.02], PlotLabel -> "FFT"];

s2c = dctFourier[s2];

ListPlot[s2c,PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"DCT"]

ListPlot[s2c[[Range[1,50]]],PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"DCT (Close-up)"]

Importing Sampled Data

Mathematica can import common file formats. AIFF (or .aif) stands for "audio interchange file format" is one of the most widely used format for storing and transmitting sampled sound. (Another format, that uses lossy data compression, is MP3 (or MPEG-1 Audio Layer-3). Mathematica can not import an MP3 file. MP3 generally reduces data size by a factor of 12 while retaining the perception of high fidelity.)

Set Directory Path

In order to import a file, the path to the directory/folder containing the data must be specified.

SetDirectory["~/Desktop/"]

Audio File Import

In order to import a file, the path to the directory/folder containing the data must be specified.

FileNames["*.aif"]

s = Import["Do_Not_Go_Gentle_64kb.aif"]

Examine Sound Object

rate = Part[s,1,2](* Sample Rate: Samples Per Second *)

1/rate // EngineeringForm(* time per sample *)

Shallow[Part[s,1,1,1]] (* Data List *)

sData = Drop[Part[s,1,1,1],1]; (* Make exactly 2^16 = 64 kS *)

nData = Length[sData]

(rate/nData) // N (* ΔFrequency Interval Hz *)

{Max[sData], Min[sData]}

ListPlot[Take[sData,{1,nData,10}],PlotRange->All,PlotLabel->"Sampled Sound Data"]

Comparing DCT and FFT

scDylan=dctFourier[sData];
ListPlot[scDylan,PlotRange->All,PlotLabel->"DCT of Dylan Thomas"]

sfDylan=Fourier[sData];
ListPlot[Abs[sfDylan],PlotRange->All,PlotLabel->"Abs[FFT] of Dylan Thomas"]

Energy DCT

eTotal = scDylan . scDylan

FoldList[Plus, 0.0, {a, b, c, d, e, f}]

eDCT = FoldList[Plus, 0, scDylan scDylan]/eTotal;

ListPlot[eDCT,PlotRange->All,PlotLabel->"Energy Fraction Integrals"]

ListPlot[eDCT[[Range[1,4000]]],PlotRange->{0,1},PlotLabel->"Energy Fraction Integrals"]

Energy Fourier

nDylan = Length[sfDylan]

sfDylanR = sfDylan[[Range[nDylan/2+1, nDylan]]] ~Join~ sfDylan[[Range[1,nDylan/2]]];

ListPlot[Abs[sfDylanR],PlotRange->All,PlotLabel->"Abs[FFT] of Dylan Thomas"]

eTotal = sfDylan . Conjugate[sfDylan] // Chop

eTotal = sfDylanR . Conjugate[sfDylanR] // Chop

eFFT = Chop[FoldList[Plus, 0, sfDylanR Conjugate[sfDylanR]]]/eTotal;

ListPlot[eFFT,PlotRange->All,PlotLabel->"Energy Integrals"]

ListPlot[eFFT[[Range[(nDylan/2)-2000,(nDylan/2)+2000]]],PlotRange->{0,1},PlotLabel->"Energy Integrals"]

Reconstructed (Filtered) Dylan

cFilter = Table[If[ 500 ≤ i ≤ 4000, 1, 0], {i,1, Length[scDylan]}];(* 300 Hz to 2.5 kHz *)

dylanC = idctFourier[cFilter scDylan];

ListPlay[dylanC,SampleRate -> rate]

ListPlay[sData,SampleRate -> rate]

2D DCT Modes

The JPEG standard is based on the DCT of 8 x 8 blocks of graphic data. The modes of an 8 x 8 block can be calculated easily.

Array[f, {8,8}]

Outer[Times, {m1, m2, m3, m4, m5}, {n1, n2, n3, n4, n5}]  // MatrixForm

dctTensor = Array[(Outer[Times, dctMatrix[8][[#1]], dctMatrix[8][[#2]]])&, {8, 8}];

Dimensions[dctTensor]

dctTensor // Shallow

ListDensityPlot[dctTensor[[2,2]],Mesh->False,PlotLabel->"DCT Mode [2,2]"]

Show[GraphicsGrid[Map[ListDensityPlot[#1,Mesh->False,Frame->False,DisplayFunction->Identity]&,dctTensor,{2}]],PlotLabel->"DCT Modes for an 8x8 Image"]

A DCT 2d Procedure

dct2D[array_?MatrixQ] := Transpose[dctFourier /@ Transpose[dctFourier /@ array]];

idct2D[array_?MatrixQ] := Transpose[idctFourier /@ Transpose[idctFourier /@ array]];

iSize=(2^6);
image1=Table[Sin[((π i 3)/iSize)] Sin[((π j 2)/iSize)],{i,1,iSize},{j,1,iSize}];
ListDensityPlot[image1,Mesh->False,PlotLabel->"Test Image (64x64)"]

image1DCT=dct2D[image1];

image1DCT // Dimensions

ListDensityPlot[image1DCT,Mesh->False,PlotRange->All,PlotLabel->"DCT of 64x64 Test Sample"]

ListPlot3D[image1DCT,PlotRange->All]

Blocked DCT 2D

JPEG encoding is performed on blocks of 8 x 8 data. In this section, a general purpose algorithm is defined that partitions a graphic into 8 x 8 blocks and then performs a DCT on the partitions. (This is based on the previously referenced article by Andrew Watson.)

Partition

Clear[a];

alist = Table[a[i], {i, 1, 36}]

Partition[alist, 8]

Blocked DCT Definitions

bDCTfactors = N[{1} ~ Join~ Table[Sqrt[2] Exp[-I π i/(2 8)], {i, 1, 8-1}]];
bIDCTfactors = Conjugate[bDCTfactors];
bDCT[list_] := Re[ bDCTfactors InverseFourier[ N[list[[{1,3,5,7,8,6,4,2}]]]]];
bIDCT[list_] := Re[Fourier[bIDCTfactors list]][[{1,8,2,7,3,6,4,5}]];
bDCT[array_?MatrixQ] := Transpose[bDCT /@ Transpose[bDCT /@ array]];
bIDCT[array_?MatrixQ] := Transpose[bIDCT /@ Transpose[bIDCT /@ array]];
bDCT[list_?(Length[#] > 8&)] := Join @@ (bDCT /@ Partition[list,8]);
bIDCT[list_?(Length[#] > 8&)] := Join @@ (bIDCT /@ Partition[list,8]);

?bDCT

?bIDCT

bimage1 = bDCT[image1];

ListDensityPlot[bimage1,Mesh->False,PlotRange->All,PlotLabel->"Block DCT of (8x8)x(8x8) Image"]

MatrixPlot[bimage1,Mesh->False,PlotRange->All,ColorFunction -> GrayLevel,ColorFunctionScaling -> False,PlotLabel->"Block DCT of (8x8)x(8x8) Image"]

Dylan Block DCT

Get Graphic Data

SetDirectory["~/Desktop/"]

FileNames["*.tif"]

g1 = Import["D_Thomas.tif"];

g1Data = ImageData[g1,"Byte"];

Dimensions[g1Data]

ListDensityPlot[g1Data,Mesh->False,AspectRatio->Automatic]

GraphicsRow[{MatrixPlot[g1Data, ColorFunction -> GrayLevel],g1Data = Reverse[g1Data,1];
MatrixPlot[g1Data, ColorFunction -> GrayLevel]}]

ListDensityPlot[g1Data,Mesh->False,AspectRatio->Automatic]

DCT

bg1 = bDCT[g1Data];
Dimensions[bg1]

ListDensityPlot[bg1,Mesh->False,AspectRatio->Automatic,PlotRange->All]

ListDensityPlot[bg1,Mesh->False,AspectRatio->Automatic]

big1Data = bIDCT[bg1];

g1Data - big1Data // Chop // Abs // Max

JPEG Quantization

In the JPEG image compression standard,each DCT coefficient is quantized using a weight that depends on the frequencies for that coefficient.The coefficients in each 8 x 8 block are divided by a corresponding entry of an 8 x 8 quantization matrix,and the result is rounded to the nearest integer.

In general,higher spatial frequencies are less visible to the human eye than low frequencies.Therefore,the quantization factors are usually chosen to be larger for the higher frequencies. (A mode is "quantized" by dividing by the quantization factor. A large number reduces the relative weight of the component.) The following quantization matrix is widely used for monochrome images and for the luminance component of a color image. According to Watson, the "de facto" quantization matrix for JPEG encoding is given by the following.

qLum =
{{16, 11, 10, 16, 24, 40, 51, 61},
{12, 12, 14, 19, 26, 58, 60, 55},
{14, 13, 16, 24, 40, 57, 69, 56},
{14, 17, 22, 29, 51, 87, 80, 62},
{18, 22, 37, 56, 68,109,103, 77},
{24, 35, 55, 64, 81,104,113, 92},
{49, 64, 78, 87,103,121,120,101},
{72, 92, 95, 98,112,100,103, 99}};

ListDensityPlot[qLum,Mesh->False,PlotRange->All,PlotLabel->"JPEG 8x8 Quantization Factors"]

Quantized Functions

BlockImage[image_, blocksize_:{8, 8}] :=Partition[image, blocksize] /;And @@ IntegerQ /@ (Dimensions[image]/blocksize)

UnBlockImage[blocks_] := Partition[Flatten[Transpose[blocks, {1, 3, 2}]],{Times @@ Dimensions[blocks][[{2, 4}]]}]

(tImage = Table[i + 8(j-1), {j, 4}, {i, 6}]) // MatrixForm

BlockImage[tImage, {2, 3}] // MatrixForm

UnBlockImage[BlockImage[tImage, {2, 3}]] // MatrixForm

bQuantize[image_, qMatrix_] := Map[(#/qMatrix)&, BlockImage[image, Dimensions[qMatrix]],{2}] //UnBlockImage // Round;

bUnQuantize[image_, qMatrix_] := Map[(# qMatrix)&, BlockImage[image, Dimensions[qMatrix]],{2}] //UnBlockImage

Examples

{Dimensions[g1Data], Dimensions[g1Data]/8}

bg1=bDCT[g1Data];
ListDensityPlot[bg1,Mesh->False,AspectRatio->Automatic,Frame->False]

bqg1=bQuantize[bg1,qLum];
ListDensityPlot[bqg1,Mesh->False,AspectRatio->Automatic,Frame->False]

buqg1=bUnQuantize[bqg1,qLum];
ListDensityPlot[buqg1,Mesh->False,AspectRatio->Automatic,Frame->False]

lp1=ListPlot[bg1[[Range[1,80],9]],Joined->True,PlotRange->{-10,256},PlotStyle->{Thickness[0.01`],RGBColor[1,0,0]},DisplayFunction->Identity];
lp2=ListPlot[buqg1[[Range[1,80],9]],Joined->True,PlotRange->{-10,256},PlotStyle->RGBColor[0,0,1],DisplayFunction->Identity];
Show[lp1,lp2,DisplayFunction->$DisplayFunction,PlotLabel->"Comparison between (Un)Quantized Data"]

lp1=ListPlot[bg1[[Range[1,80],9]],Joined->True,PlotRange->{-10,25},PlotStyle->{Thickness[0.01`],RGBColor[1,0,0]},DisplayFunction->Identity];
lp2=ListPlot[buqg1[[Range[1,80],9]],Joined->True,PlotRange->{-10,25},PlotStyle->RGBColor[0,0,1],DisplayFunction->Identity];
Show[lp1,lp2,DisplayFunction->$DisplayFunction,PlotLabel->"Comparison between (Un)Quantized Data"]

Data Compression

Count[Flatten[bg1], _?(# == 0 &)]

Count[Flatten[buqg1], _?(# == 0 &)]

Histogram[Flatten[bg1],PlotLabel->"DCT Frequencies"]

Histogram[Flatten[buqg1],PlotLabel->"Quantized Frequencies"]

Default Compression

q1Data = g1Data // bDCT // bQuantize[#,qLum]& // bUnQuantize[#, qLum]& // bIDCT;

Show[GraphicsRow[(ListDensityPlot[#1,PlotRange->{0,255},Mesh->False,AspectRatio->Automatic,DisplayFunction->Identity,Frame->False]&)/@{g1Data,q1Data}]]

Low Quality Compression

q1Data8 = g1Data // bDCT // bQuantize[#,8qLum]& // bUnQuantize[#, 8qLum]& // bIDCT;

Show[GraphicsRow[(ListDensityPlot[#1,PlotRange->{0,255},Mesh->False,AspectRatio->Automatic,DisplayFunction->Identity,Frame->False]&)/@{g1Data,q1Data8}]]

High Quality Compression

q1Data025 = g1Data // bDCT // bQuantize[#,0.25qLum]& // bUnQuantize[#, 0.25qLum]& // bIDCT;

Show[GraphicsRow[(ListDensityPlot[#1,PlotRange->{0,255},Mesh->False,AspectRatio->Automatic,DisplayFunction->Identity,Frame->False]&)/@{g1Data,q1Data025}]]

Summary

This notebook defines and illustrates the use of the discrete Cosine transform (DCT) for digital signal processing. The DCT transforms a sequence of Real numbers into another sequence of Real numbers. Since it is related to the Fourier transform (DFT), relatively fast algorithms exist that have been implemented both in hardware and software.

When the DCT is applied to two-dimensional images, digital filters and compression technqiues can be applied. The JPEG algorithm is applied to 8 x 8 blocks using a perception-based quantization. The quality of the image and the degree of compression can be adjusted with a single parameter.