The Discrete Cosine Transform and Chebyshev Expansions

Cosine series and the discrete Cosine transform are close relatives of Fourier series and the discrete Fourier transform. A Fourier series in triginometric polynomials is an expansion

$$f(\theta) = \sum_{k=0}^\infty c_k \cos k \theta + \sum_{k=1}^\infty s_k \sin k \theta,$$

where

$$ \begin{align*} c_0 &= {1 \over 2 \pi} \int_0^{2 \pi} f(\theta) \cos k \theta d \theta \cr c_k &= {1 \over \pi} \int_0^{2 \pi} f(\theta) \cos k \theta d \theta \cr s_k &= {1 \over \pi} \int_0^{2 \pi} f(\theta) \sin k \theta d \theta \cr \end{align*} $$

A Cosine series is an expansion like

$$f(\theta) = \sum_{k=0}^\infty c_k \cos k \theta$$

That is, it is a Fourier series with only cosine terms, no sine terms.

Cosine series for even functions

Even functions satisfy $f(\theta) = f(-\theta)$. Cosine expansions capture evenness, because each Cosine term is even: $\cos k \theta = \cos(-k\theta)$.

Exercise 1(a) Why is it true that, if $f$ is even, we have

$$s_k = 0?$$

Verify this statement for $k=1,2,3$ using quadgk for the even function $f(\theta) = e^{\cos \theta}$.

In [1]:
f=θ->exp(cos(θ))

quadgk(θ->f(θ)*sin(θ),0,2π)[1],quadgk(θ->f(θ)*sin(2θ),0,2π)[1],quadgk(θ->f(θ)*sin(3θ),0,2π)[1]
Out[1]:
(-3.80064384882164e-16,3.674734011551995e-15,2.651001459313375e-15)

Exercise 1(b) When $f$ is even, the integrals can be reduced to have a period. Explain why its true that

\begin{align*} c_0 &= {1 \over \pi} \int_0^{\pi} f(\theta) \cos k \theta d \theta \cr c_k &= {2\over \pi} \int_0^{\pi} f(\theta) \cos k \theta d \theta \cr \end{align*}

Verify this statement using quadgk for the even function $f(\theta) = e^{\cos \theta}$ for $k=0,1,2,3$.

In [2]:
quadgk(θ->f(θ),0,2π)[1]/(2π)-quadgk(θ->f(θ),0,π)[1]/(π)
k=1;quadgk(θ->f(θ)*cos(k*θ),0,2π)[1]/(π)-2quadgk(θ->f(θ)*cos(k*θ),0,π)[1]/(π)
k=2;quadgk(θ->f(θ)*cos(k*θ),0,2π)[1]/(π)-2quadgk(θ->f(θ)*cos(k*θ),0,π)[1]/(π)
k=3;quadgk(θ->f(θ)*cos(k*θ),0,2π)[1]/(π)-2quadgk(θ->f(θ)*cos(k*θ),0,π)[1]/(π)
Out[2]:
-9.159339953157541e-16

Exercise 1(c) Use quadgk to approximate $c_k$ for $f(\theta) = e^{\cos \theta}$. Verify

$$f(\theta) \approx \sum_{k=0}^9 c_k \cos k \theta$$
In [3]:
c=Array(Float64,10)
k=0;c[1]=quadgk(θ->f(θ)*cos(k*θ),0,π)[1]/(π)
for k=1:length(c)-1
    c[k+1]=2quadgk(θ->f(θ)*cos(k*θ),0,π)[1]/(π)
end


t=0.1
ret=0.

for k=0:length(c)-1
    ret+=c[k+1]*cos(k*t)
end
ret-f(t)
Out[3]:
-3.091917832875879e-10

Non-periodic functions and Chebyshev expansion

Consider a non-periodic function $f(x)$ on $[-1,1]$. Then $f(\cos \theta)$ is a periodic, even function on $[0,2\pi]$! Cosine expansion of $f(\cos \theta)$ combined with the change of variables $\theta = \hbox{acos}\, x$ leads to the Chebyshev expansion of a non-periodic function. Remarkably, this expansion is a polynomial expansion.

Exercise 2(a) Consider $f(x) = e^x$. Show that

$$f(x) \approx \sum_{k=0}^9 c_k T_k(x)$$

where $c_k$ are the cosine coefficients of $f(\cos \theta)$ and $T_k(x) \triangleq \cos k \, {\rm acos}\, x$.

In [ ]:

Exercise 2(b) Advanced Use the fact that

$$\int_0^\pi \cos k \theta \sin \theta d \theta = \begin{cases} 2/(1-k^2) & \hbox{$k$ even} \cr 0 & \hbox{otherwise} \end{cases}$$

to approximate

$$\int_{-1}^1 f(x) dx \approx \int_{-1}^1 \sum_{k=0}^{n-1} c_k T_k(x) dx$$

This is known as Clenshaw–Curtis quadrature.

In [ ]:

Exercise 2(c) Advanced Show that $T_k(x)$ are polynomials.

Exercise 2(d) Advanced Show that

$$\int_{-1}^1 {T_k(x) T_j(y) \over \sqrt{1-x^2}} dx = 0 $$

if $k \neq j$.

In [ ]:

The DCT

The Trapezium rule over $[0,\pi]$ also leads to a discrete cosine transform. We can use the trap routine from lectures:

In [1]:
function trap(f,a,b,n)
    h=(b-a)/n
    x=linspace(a,b,n+1)

    v=f(x)
    h/2*v[1]+sum(v[2:end-1]*h)+h/2*v[end]
end
Out[1]:
trap (generic function with 1 method)

Exercise 3(a) Use trap to complete the function mydct that returns approximations $[c_0^n,\ldots,c_{n-1}^n]$ to $c_k$. (Hint: make sure you use n-1 when calling trap.) Show experimentally that it is exact for $f(\theta) = \cos \omega \theta$ when $\omega > n$.

In [17]:
function mydct(f::Function,n)
    c=zeros(n)
    c[1]=trap(f,0,π,n)/π
    for k=1:n-1
        c[k+1]=trap(θ->cos(k*θ).*f(θ),0,π,n)*2/π
    end
    c
end
Out[17]:
mydct (generic function with 1 method)
In [30]:
mydct(θ->cos(13θ),10)
Out[30]:
10-element Array{Float64,1}:
  3.53395e-16
 -4.59413e-16
  0.0        
  5.30092e-16
 -6.7145e-16 
 -1.06018e-16
  1.76697e-17
  1.0        
  1.76697e-17
  1.41358e-16

The following command is a fast DCT-I.

In [6]:
dctI(v) = FFTW.r2r(v, FFTW.REDFT00)
Out[6]:
dctI (generic function with 1 method)

Exercise 3(b) Advanced Figure out how to calculate $c_k^n$ using dctI and $f$ evaluated at θ=linspace(0,π,n). (Hint: try it with $\cos k \theta$ for different choices of $k$.)

In [7]:
n=11;θ=linspace(0,π,n)
v=exp(cos(θ-0.1))
cfs=dctI(v)/(n-1)
cfs[1]/=2
cfs[end]/=2

cfs
Out[7]:
11-element Array{Float64,1}:
  1.33991   
  1.16995   
  0.226631  
  0.0151554 
 -0.0132564 
 -0.00865769
 -0.00802881
 -0.00500644
 -0.00554588
 -0.00393199
 -0.00248131
In [ ]: