Issue
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.
Answered By - Nick ODell
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.