1

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.pik(2n[0]+1))/(2N)) + 
    y[1] * np.cos((np.pik(2n[1]+1))/(2N)) + 
    y[2] * np.cos((np.pik(2n[2]+1))/(2N)))

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)

Antoni Parellada
  • 453
  • 4
  • 14

1 Answers1

1

I'd take bets that the 2D-DCT is just separable into identical row-wise and column-wise DCTs. Anything else would be surprising.

Anyways, you know the 2D-DCT is a reversible linear operation: It maintains dimensionality, and has only a trivial null space. Therefore, you can analyze what it does directly from applying it to an arbitrarily chosen basis of your matrix space.

Marcus Müller
  • 30,525
  • 4
  • 34
  • 58
  • 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? – Antoni Parellada May 05 '21 at 21:04
  • You want to learn what a linear operator (the 2D DFT) does. So, do what you'd do to figure that out: Test it on a set of vectors that form a base. – Marcus Müller May 05 '21 at 22:25