[PPT] - Algorithms of Scientific Computing Discrete Cosine Transform (DCT) PowerPoint Presentation

SLIDE 1

Technische Universit¨ at M¨ unchen

Algorithms of Scientific Computing

Discrete Cosine Transform (DCT)

Michael Bader

Summer Term 2012

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 1

SLIDE 2

Technische Universit¨ at M¨ unchen

DFT and Symmetry

INPUT TRANSFORM real symmetry fn ∈ R → Real DFT (RDFT) even symmetry fn = f−n → Discrete Cosine Transform (DCT)

dd symmetry

fn = −f−n → Discrete Sine Transform (DST) “QUARTER-WAVE” INPUT TRANSFORM even symmetry fn = f−n−1 → QW-DCT

dd symmetry

fn = −f−n−1 → QW-DST

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 2

SLIDE 3

Technische Universit¨ at M¨ unchen

Application Example: Compression of Image Data (JPEG)

Compression steps of the JPEG method

1. Conversion into a suitable colour model (YCbCr, e.g.),

separation of brightness and colour information

2. Downsampling (in particular of the colour components)
3. blockwise “quarter-wave discrete cosine transform”

(blocks of size 8 × 8)

4. Quantification of the coefficients (→ reduce information)
5. run-length encoding, Huffman/arithmetic coding

(loss-free compression of the quantified coefficients) Example: jpeg on matlab central (see link on webpage)

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 3

SLIDE 4

Technische Universit¨ at M¨ unchen

Discrete Fourier Transform (DFT)

Definition: For a vector of N complex numbers (f0, . . . , fN−1), the discrete Fourier transform is given by the vector (F0, . . . , FN−1), where Fk = 1 N

N−1

n=0

fne−i2πnk/N. Interpretation:

as trigonometric interpolation/approxmiation
as approximation of the coefficients of the Fourier series

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 4

SLIDE 5

Technische Universit¨ at M¨ unchen

Fourier Coefficients and Numerical Quadrature

For a 2π-periodic function f, the corresponding Fourier series is defined as f(x) ∼

∞

k=−∞

ckeikx, mit ck = 1 2π

2π

f(x)e−ikx dx

The ck are called (continuous) Fourier coefficients. If f is piecewisely smooth, Fourier series converges pointwisely (i.e. for each x) towards

1 2(f(x+) + f(x−)),

i.e. in particular towards f(x), if f is coninuously differentiable at x.

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 5

SLIDE 6

Technische Universit¨ at M¨ unchen

Approximate Computation of ck

The continuous Fourier coefficients are given as ck = 1 2π

2π

f(x)e−ikx dx

Approaches to compute ck approximately:

compute ck only for ±k = 0, . . . , K; then: f(x) ≈

K

k=−K

ckeikx

compute integral

2π

f(x)e−ikx dx numerically

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 6

SLIDE 7

Technische Universit¨ at M¨ unchen

Computation of ck via Trapezoidal Sum

Trapezoidal sum : for equidistant xn := 2πn

N : 2π

g(x) dx ≈ TN{g} := 2π

N

1

2g(x0) +

N−1

n=1

g(xn) + 1 2g(xN)

Use g(x) := f(x)e−ikx and fn := f(xn); hence:

ck ≈ 1 2π TN{f(x)e−ikx} = 1 N

1

2f0e0 +

N−1

n=1

fne−i2πnk/N + 1 2fNe−i2πNk/N

=

1 N

f0

2 +

N−1

n=1

fne−i2πnk/N + fN 2

Michael Bader: Algorithms of Scientific Computing

Discrete Cosine Transform (DCT), Summer Term 2012 7

SLIDE 8

Technische Universit¨ at M¨ unchen

Computation of ck via Trapezoidal Sum (2)

If f0 = fN (periodic data), we obtain ck ≈ Fk = 1 N

N−1

n=0

fne−i2πnk/N ⇒ Fk are approximations of ck ⇒ approximate computation leads to solution of the interpolation problem ⇒ approximation error is of order O(N−2) For f0 = fN, or for “discontinuities”, we get a recommendation: Average Values at Endpoints and Discontinuities (AVED)

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 8

SLIDE 9

Technische Universit¨ at M¨ unchen

Computation of ck via Midpoint Rule

Midpoint rule: evaluate g(x) at midpoints xn:

2π

g(x) dx ≈ 2π

N

N−1

n=0

g(xn) with xn := 2π

n + 1

2

N

. With g(x) := f(x)e−ikx and fn := f(xn), we obtain: ck ≈ Fk := 1 N

N−1

n=0

fne−i2π(n+ 1

2)k/N

“Quarter-Wave Discrete Fourier Transform”

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 9

SLIDE 10

Technische Universit¨ at M¨ unchen

Quarter-Wave Discrete Fourier Transform

new variant of DFT:
Fk := 1

N

N−1

n=0

fne−i2π(n+ 1

2)k/N

fn :=

N−1

k=0
Fkei2π(n+ 1

2)k/N

Comparison with coefficients Fk of the “usual” DFT:

Fk = Fkeiπk/N = Fkωk/2

N

Supporting points compared to “usual” DFT shifted by a “quarter

wave length” (midpoints of intervals).

Derivation via midpoint rule motivates usage for piecewise

constant data ⇒ Transformation of image data

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 10

SLIDE 11

Technische Universit¨ at M¨ unchen

Quarter-Wave DFT on Symmetric Data

Given 2N real-valued input data f0, . . . , f2N−1 with symmetry f2N−n−1 = fn Inserting the symmetric data in Quarter-Wave DFT results in

Fk

= 1 2N

2N−1

n=0

fnω

−k(n+ 1

2)

2N

= 1 2N

N−1

n=0

fnω

−k(n+ 1

2)

2N

+ 1 2N

N−1

n=0

f2N−n−1ω

−k(2N−n−1+ 1

2)

2N

= 1 2N

N−1

n=0

fn

ω

−k(n+ 1

2)

2N

+ ω

−k(−n− 1

2)

2N

= 1

N

N−1

n=0

fn cos

πk
n + 1

2

N
.

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 11

SLIDE 12

Technische Universit¨ at M¨ unchen

Quarter-Wave DFT on Symmetric Data (2)

Quarter-Wave DFT of symmetric data results in real-valued coefficients:

Fk = 1

N

N−1

n=0

fn cos

πk
n + 1

2

N
for k = 0, . . . , 2N − 1

Additional symmetry:

F2N−k

= 1 N

N−1

n=0

fn cos

π(2N − k)
n + 1

2

N
=

1 N

N−1

n=0

fn cos

2πn + π − πk
n + 1

2

N
= −

Fk ⇒ again: only N independent coefficients

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 12

SLIDE 13

Technische Universit¨ at M¨ unchen

Quarter-Wave Even Discrete Cosine Transform

backward transform: fn :=

2N−1

k=0
Fkei2π(n+ 1

2)k/2N

F2N−k=−

Fk

− → fn = F0+2

N−1

k=1
Fk cos
πk
n + 1

2

N
Definition of the quarter-wave even DCT:
Fk = 1

N

N−1

n=0

fn cos

πk
n + 1

2

N
fn =

F0+2

N−1

k=1
Fk cos
πk
n + 1

2

N
N real-valued data

← → N real-valued coefficients (no symmetry any more in data/coefficients!)

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 13

SLIDE 14

Technische Universit¨ at M¨ unchen

2D Cosine Transform

Definition of the 2D-DCT:

Fkl

= 1 N · M

N−1

n=0

M−1

m=0

fnm cos

πk
n + 1

2

N
cos
πl
m + 1

2

M
fnm

= 4

N−1

k=0

′ M−1

l=0

′

Fkl cos

πk
n + 1

2

N
cos
πl
m + 1

2

M
shortened notation:

N−1

k=0

′ xk := x0 2 + N−1

k=1

xk Application: blockwise 2D-DCT in JPEG/MPEG compression

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 14

SLIDE 15

Technische Universit¨ at M¨ unchen

Reduction of the 2D-FCT to 1D-FCTs

In the 2D cosine transform, we can rearrange:

Fkl

= 1 N2

N−1

n=0

N−1

m=0

fnm cos

πk
n + 1

2

N
cos
πl
m + 1

2

N
=

1 N

N−1

n=0
1

N

N−1

m=0

fnm cos

πl
m + 1

2

N
:=

Fnl cos

πk
n + 1

2

N
.
For each n,

Fnl are computed via N 1D transforms

we may first 1D-transform all rows and then all columns to get

the 2D-transform

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 15

SLIDE 16

Technische Universit¨ at M¨ unchen

Application Example: Compression of Image Data (JPEG)

Compression steps of the JPEG method

1. Conversion into a suitable colour model (YCbCr, e.g.),

separation of brightness and colour information

2. Downsampling (in particular of the colour components)
3. blockwise “quarter-wave discrete cosine transform”

(blocks of size 8 × 8)

4. Quantification of the coefficients (→ reduce information)
5. run-length encoding, Huffman/arithmetic coding

(loss-free compression of the quantified coefficients) Example: jpeg on matlab central (see link on webpage)

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 16

SLIDE 17

Technische Universit¨ at M¨ unchen

QW-DCT – Algorithm

Reduce to Real FFT: (1) for n = 0, . . . , N − 1: gn = fn g2N−n−1 = fn (2) 2N-Real-FFT: compute Gk from gn (for k = 0, . . . , N) (3) for k = 0, . . . , N − 1:

Fk = Gke−iπk/2N

Further optimisations:

substitute real-valued 2N-FFT by complex N-FFT
compact (divide-and-conquer) real FFT
compact Fast DCT −

→ paper Swarztrauber

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 17

SLIDE 18

Technische Universit¨ at M¨ unchen

Compact Fast DCT

Consider DCT: with symmetry f2N−n−1 = fn

Fk = 1

2N

2N−1

n=0

fnω

−k(n+ 1

2)

2N

− →

Fk = 1

N

N−1

n=0

fn cos

πk
n + 1

2

N
.

Split into even and odd indices: gn := f2n and hn := f2n+1 (as in FFT)

gn := f2n:

gn = f2n = f2N−2n−1 = f2(N−n)−1 = f2(N−n−1)+1 = hN−n−1

hn := f2n+1:

hn = f2n+1 = f2N−(2n+1)−1 = f2(N−n−1) = gN−n−1

thus: two real DFTs with symmetric data sets

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 18

SLIDE 19

Technische Universit¨ at M¨ unchen

Compact Fast IDCT

Consider backward transform: with symmetry F2N−k = − Fk fn :=

2N−1

k=0
Fkei2π(n+ 1

2)k/2N

− → fn = F0+2

N−1

k=1
Fk cos
πk
n + 1

2

N
Split into even and odd indices: (as in FFT)
Gk :=

F2k: leads to IDCT −Gk = − F2k = F2N−2k = F2(N−k) = GN−k

Hk :=

F2k+1: leads to “DCT with negative exponent”?? −Hk = − F2k+1 = F2N−(2k+1) = F2(N−k)−1 = F2(N−k−1)+1 = HN−k−1 − → to be continued . . .

Michael Bader: Algorithms of Scientific Computing Discrete Cosine Transform (DCT), Summer Term 2012 19