Applying DCT matrix along each axis not giving the desired result

0
votes

Using the following code, I would expect to get linalg.norm(y-z) equals zero. I follow the idea of Apply DFT matrix along each axis of 3D array in NumPy? since FFT and DCT are both separable, unitary linear transforms. But somehow this isn't the case.

import numpy as np
from scipy.fftpack import dct 

x = np.random.rand(5,5)

D = dct(np.eye(5), norm='ortho')
y = np.dot(D,np.dot(x, D.T))
z = dct(dct(x, axis = 0 , norm = 'ortho'), axis = 1 , norm = 'ortho')
1
pythonnumpyfftdct
@Ehsan Why did you change that to DFT? Isn't the question clearly about DCT? – Mateen Ulhaq
@MateenUlhaq My bad. Misunderstood it. Thank you for the correction. – Ehsan

1 Answers

2
votes

Pick one of the following solutions.

  1. Redefine D by specifying axis=0:

    D = dct(np.eye(n), axis=0, norm="ortho")

  2. Redefine D by specifying .T and using the default axis=-1 (doesn't really generalize to higher dimensions...):

    D = dct(np.eye(n), norm="ortho").T

  3. Use D.T @ x to represent a DCT along axis=0 of x:

    y = D.T @ x @ D

The order doesn't matter in the case of DFT since the matrix is symmetric (D == D.T) in addition to being unitary (D.conj().T @ D == 1). But the DCT matrix is not symmetric, so you have to be careful of what order you use.

Consider the definition for DCT-II:

enter image description here

When you construct the operator D to imply that D @ x takes the DCT-II along the rows (axis=0), then D must be defined so that the covariance and contravariance of the transformation are properly handled.

Full example:

import numpy as np
from scipy.fftpack import dct 

n = 5
x = np.random.rand(n, n)
D = dct(np.eye(n), axis=0, norm="ortho")

y = D @ x @ D.T

z = x
z = dct(z, axis=0, norm="ortho")
z = dct(z, axis=1, norm="ortho")

>>> np.linalg.norm(z - y)
6.20161216470283e-16