# How can extract the cosine transform formula used for 2D by scipy.fft.dct

For the 1D cosine transform the documentation is clear in here, and I can reproduce it easily:

The formula is:

$$y_k= 2 \sum_{n=0}^{N-1}x_n \cos \left( \frac{\pi k(2n + 1)}{2N}\right)$$

and here it is reproduced manually, for example for the third harmonic:

import numpy as np
from scipy.fft import fft, dct

y = np.array([10,3,5])

# dct() call:
print(dct(y, 2))

# manual reproduction for k = 2:
N = 3
k = 2
n = np.array([0,1,2])
2 * (
y[0] * np.cos((np.pi*k*(2*n[0]+1))/(2*N)) +
y[1] * np.cos((np.pi*k*(2*n[1]+1))/(2*N)) +
y[2] * np.cos((np.pi*k*(2*n[2]+1))/(2*N)))


in both cases I get a 9.

But there is no documentation for the 2D DCT, and my attempts at hacking the formula with a toy matrix are not working out:

Compare:

z = np.array([[ 2,  3       ],
[ 10, 15]])
dct(z, axis=0)  # dct() call


to for instance:

N = 2
M = 2
k = 0
l = 0
n = np.array([0,1])
m = np.array([0,1])
M*N * (
z[0,0] * np.cos((np.pi*k*(2*n[0]+1))/(2*N)) * np.cos((np.pi*l*(2*m[0]+1))/(2*M)) +
z[0,1] * np.cos((np.pi*k*(2*n[0]+1))/(2*N)) * np.cos((np.pi*l*(2*m[1]+1))/(2*M)) +
z[1,0] * np.cos((np.pi*k*(2*n[1]+1))/(2*N)) * np.cos((np.pi*l*(2*m[0]+1))/(2*M)) +
z[1,1] * np.cos((np.pi*k*(2*n[1]+1))/(2*N)) * np.cos((np.pi*l*(2*m[1]+1))/(2*M))
)


for the first coefficient.

Can anyone help me reconcile the output of dct() with my attempted manual calculation?

I guess the formula is not...

$$y_{k,l}= 4 \sum_{n=0}^{N-1}\sum_{m=0}^{M-1}x_n y_m \cos \left( \frac{\pi k(2n + 1)}{2N}\right)\cos \left( \frac{\pi l(2m + 1)}{2M}\right)$$

but it would be really easy to correct if I could get the same output manually for one of the coefficients on the example matrix above.

The methods applied here make explicit reference to this author, but they don't help me understand what the function is doing.

PS The formula turns out to be spot on if it is called twice as in:

dct(dct(z.T, axis=0).T, axis=0)

• Don't you just apply the DCT in one dimension, then again in the other? – endolith May 5 at 19:36

• Do you mean in your last paragraph that applying the dct() call twice to the identity matrix as in dct(dct(np.eye(4).T, axis=0).T, axis=0) would make sense? What could I learn from it? – Math tourist May 5 at 21:04