Poor man’s discrete cosine transform

The discrete cosine transform can be understood as a trigonometric interpolation problem. Given some data, like {y = (3,1,4,1,5,9)}, I want to find the coefficients {A_0,\dots,A_5} such that

\displaystyle \sum_{k=0}^5 A_k\cos kx_j = y_j,\quad j= 0,\dots,5

with {x_j=\pi j/6}. (This is different from the unitary scaling preferred in image processing.) Matlab’s signal processing toolbox, as well as “signal” package of Octave-Forge both have dcp command. But the former is expensive and the latter is not entirely painless to install. Let’s work with the standard tool fft, readily available in Matlab, Scilab, and Octave.

Screenshot 2015-12-13 at 2.37.29 AM
Points to interpolate

The full Fourier series reduces to cosines when the function is even. Therefore, we need to extend {y} by an even reflection. One might imagine such an extension as

\displaystyle (3,1,4,1,5,9,9,5,1,4,1,3)

but this is incorrect. Thinking of the six {y} values as the samples of a function at {0,\pi/5,2\pi/5,\dots,\pi} helps: the even reflection across {\pi} should not duplicate the boundary value at {\pi}. Also, {0} gets reflected to {2\pi}, which is the same as {0} due to periodicity, so this value also should not be repeated. The appropriate reflection is

\displaystyle (3,1,4,1,5,9,5,1,4,1)

In Matlab terms,

y = [3 1 4 1 5 9]
z = [y, y(end-1:-1:2)]
w = fft(y)

The output is real, as desired. It also has the same kind of symmetry as the input.

\displaystyle (34, -10.618, 7.618, -8.382, 5.382, 8, 5.382, -8.382, 7.618, -10.618)

The term {34} is simply the sum of the {z}-values, but the constant term in Fourier expansion is the mean value. So, division by the data size ({10}) is needed to get the coefficients. Now that we have them, it seems that the interpolating function should be

\displaystyle 3.4 -1.0618\cos x+0.7618 \cos 2x -\dots -1.0618\cos 9x

but this is incorrect. Although this function does pass through the given points,

Screenshot 2015-12-13 at 2.38.24 AM
Interpolation, 1st attempt

it wiggles too much due to high frequency terms such as {\cos 9x}. The aliasing comes into play: at the multiples of {\pi/5}, the function {\cos 9x} is the same as {\cos x}.

Screenshot 2015-12-13 at 2.39.32 AM
Aliasing: at the sample points marked in red, the functions agree

Similarly, {\cos 8x} should be replaced by {\cos 2x} in the interpolating function, etc. This results in a much smoother interpolant, which is a linear combination of cosines up to {\cos 5x}.

Screenshot 2015-12-13 at 2.48.08 AM.png
Correct cosine expansion in green

The final answer is

\displaystyle y= 3.4 -2.1236\cos x + 1.5236\cos 2x -1.6764\cos 3x +1.0764\cos 4x + 8\cos 5x

Here is the complete Matlab code for this process of reflection, transforming and subsequent “folding” of the coefficients.

y = [3 1 4 1 5 9]
z = [y, y(end-1:-1:2)]
w = real(fft(z))/length(z)
a = [w(1), 2*w(2:end/2), w(end/2+1)]

Taking real part of FFT makes sure that nothing imaginary slips in through the errors associated with numeric computation.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s