pythonnumpyscipydct

How to "de-orthogonalize" type 2 DCT result


I have been given the task of rewriting some python code into a different language. A piece of that codes makes use of the standard 2-D DCT in scipy. I honestly don't have much clue about different types of DCT beside the very high-level idea. So I tried to implement the 2-D DCT based on the simple definition here.

Now the result of doing that does not match the default scipy dctn result. Only if you specify norm="ortho" as a parameter does the result match. Is there an easy operation I can do either before or afterwards to make the results match what dctn gives?


import numpy as np
import scipy.fftpack

def dct_mat(data):
    (n, m) = data.shape
    assert(m == n)

    dv = 1.0 / np.sqrt(float(n))
    c1 = np.sqrt(2.0 / float(n))

    mat = np.zeros((n, n))

    for i in range(n):
        mat[0][i] = dv

    for r in range(1, n):
        for c in range(n):
            mat[r][c] = c1 * np.cos((np.pi * (2.0 * float(c) + 1.0) * float(r)) / (2.0 * float(n)))

    return mat

test_data = np.array([
    [13.4, 13.9, 0.1],
    [20.1, -20.3, 88.2],
    [0, 4.2, -53.5]
])

scipy = scipy.fftpack.dctn(test_data) # with norm="ortho" arg it matches

mat = dct_mat(test_data)
res = np.matmul(np.matmul(mat, test_data), mat.T)

print("Scipy result:")
print(scipy)
print("my result:")
print(res)

Solution

  • I find that it provides the same results as scipy's non-normalized method if I change it like so:

    def dct_mat(data):
        (n, m) = data.shape
        assert(m == n)
    
        dv = 2
        c1 = 2
    
        mat = np.zeros((n, n))
    
        for i in range(n):
            mat[0][i] = dv
    
        for r in range(1, n):
            for c in range(n):
                mat[r][c] = c1 * np.cos((np.pi * (2.0 * float(c) + 1.0) * float(r)) / (2.0 * float(n)))
    
        return mat
    

    Don't ask me why this works; I don't know.

    I tested this on your original test case plus a few randomly generated matrices.