# 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.

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,

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}$.

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}$.

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.