CODEWITHSUNDEEP
pythonscipyfft
Antoni Parellada
Antoni Parellada

Reputation: 4791

How can I 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:

screenshot from scipy.fft.dct documentation detailing Type II DCT

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

screenshot of a generalized 2d DCT formula, crossed out with a red line

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.

Upvotes: 2

Views: 659

Answers (1)

Andras Deak -- Слава Україні
Andras Deak -- Слава Україні

Reputation: 35080

You can't find the formula for the multidimensional mode because the function doesn't do multidimensional cosine transforms. The axis keyword should be suspicious: in NumPy, SciPy it typically determines the direction along which lower-dimensional operations should be performed.

In other words, dct(z, axis=0) is just a columnwise 1d cosine transform:

import numpy as np
from scipy.fft import dct

z = np.array([[ 2,  3],
             [ 10, 15]])
print(dct(z, axis=0))  # dct() call
print(np.array([dct(column) for column in z.T]).T)
# both outputs
# [[ 24.          36.        ]
#  [-11.3137085  -16.97056275]]

Note the two transposes on the last line: all that it did was first loop over the array to slice it according to columns, then join them again columnwise. The latter would probably be better spelled out as

res = np.stack([dct(column) for column in z.T], axis=1)

Upvotes: 1

Related Questions