Scipy's fftpack dct and idct

Question

Let say you use the dct function, then do no manipulation of the data and use the invert transform; wouldn't the inverted data be the same as the pre-transformed data? Why the floating point issue? Is it a reported issue or is it a normal behavior?

In [21]: a = [1.2, 3.4, 5.1, 2.3, 4.5]

In [22]: b = dct(a)

In [23]: b
Out[23]: array([ 33.        ,  -4.98384545,  -4.5       ,  -5.971707  ,   4.5       ])

In [24]: c = idct(b)

In [25]: c
Out[25]: array([ 12.,  34.,  51.,  23.,  45.])

Anyone has an explanation as why? Of course, a simple c*10**-1 would do the trick, but if you repeat the call of the function to use it on several dimensions, the error gets bigger:

In [37]: a = np.random.rand(3,3,3)

In [38]: d = dct(dct(dct(a).transpose(0,2,1)).transpose(2,1,0)).transpose(2,1,0).transpose(0,2,1)

In [39]: e = idct(idct(idct(d).transpose(0,2,1)).transpose(2,1,0)).transpose(2,1,0).transpose(0,2,1)

In [40]: a
Out[40]: 
array([[[ 0.48709809,  0.50624831,  0.91190972],
        [ 0.56545798,  0.85695062,  0.62484782],
        [ 0.96092354,  0.17453537,  0.17884233]],

       [[ 0.29433402,  0.08540074,  0.18574437],
        [ 0.09942075,  0.78902363,  0.62663572],
        [ 0.20372951,  0.67039551,  0.52292875]],

       [[ 0.79952289,  0.48221372,  0.43838685],
        [ 0.25559683,  0.39549153,  0.84129493],
        [ 0.69093533,  0.71522961,  0.16522915]]])

In [41]: e
Out[41]: 
array([[[ 105.21318703,  109.34963575,  196.97249887],
        [ 122.13892469,  185.10133376,  134.96712825],
        [ 207.55948396,   37.69964085,   38.62994399]],

       [[  63.57614855,   18.44656009,   40.12078466],
        [  21.47488098,  170.42910452,  135.35331646],
        [  44.00557341,  144.80543099,  112.95260949]],

       [[ 172.69694529,  104.15816275,   94.69156014],
        [  55.20891593,   85.42617016,  181.71970442],
        [ 149.2420308 ,  154.48959477,   35.68949734]]])

Here a link to the doc.

Use `norm='ortho'` to get the normalized answers. By default `dct` and `idct` do not seem to be normalizing. — Pavan Yalamanchili, Jan 14 '13 at 20:25
Thanks so much! I did not realize that I was somewhat giving the answer in my question! — macrocosme, Jan 14 '13 at 20:47

Pavan Yalamanchili · Accepted Answer · 2015-06-04T14:58:38.747

It looks like dct and idct do not normalize by default. define dct to call fftpack.dct in the following manner. Do the same for idct.

In [13]: dct = lambda x: fftpack.dct(x, norm='ortho')

In [14]: idct = lambda x: fftpack.idct(x, norm='ortho')

Once done, you will get back the original answers after performing the transforms.

In [19]: import numpy

In [20]: a = numpy.random.rand(3,3,3)

In [21]: d = dct(dct(dct(a).transpose(0,2,1)).transpose(2,1,0)).transpose(2,1,0).transpose(0,2,1)

In [22]: e = idct(idct(idct(d).transpose(0,2,1)).transpose(2,1,0)).transpose(2,1,0).transpose(0,2,1)

In [23]: a
Out[23]: 
array([[[ 0.51699637,  0.42946223,  0.89843545],
        [ 0.27853391,  0.8931508 ,  0.34319118],
        [ 0.51984431,  0.09217771,  0.78764716]],

       [[ 0.25019845,  0.92622331,  0.06111409],
        [ 0.81363641,  0.06093368,  0.13123373],
        [ 0.47268657,  0.39635091,  0.77978269]],

       [[ 0.86098829,  0.07901332,  0.82169182],
        [ 0.12560088,  0.78210188,  0.69805434],
        [ 0.33544628,  0.81540172,  0.9393219 ]]])

In [24]: e
Out[24]: 
array([[[ 0.51699637,  0.42946223,  0.89843545],
        [ 0.27853391,  0.8931508 ,  0.34319118],
        [ 0.51984431,  0.09217771,  0.78764716]],

       [[ 0.25019845,  0.92622331,  0.06111409],
        [ 0.81363641,  0.06093368,  0.13123373],
        [ 0.47268657,  0.39635091,  0.77978269]],

       [[ 0.86098829,  0.07901332,  0.82169182],
        [ 0.12560088,  0.78210188,  0.69805434],
        [ 0.33544628,  0.81540172,  0.9393219 ]]])

I am not sure why no normalization was chosen by default. But when using ortho, dct and idct each seem to normalize by a factor of 1/sqrt(2 * N) or 1/sqrt(4 * N). There may be applications where the normalization is needed for dct and not idct and vice versa.

I just realized that it does not lead to a perfect inversion: In [185]: c-a Out[185]: array([[[ 0.00000000e+00, -2.22044605e-16, -3.33066907e-16],[ -1.11022302e-16, -3.88578059e-16, -1.11022302e-16],[ -3.33066907e-16, -1.11022302e-16, -5.55111512e-17]],[[ -3.33066907e-16, -4.44089210e-16, 0.00000000e+00],[ -1.11022302e-16, -1.11022302e-16, -5.55111512e-16],[ -2.22044605e-16, -6.66133815e-16, -3.33066907e-16]],[[ -3.33066907e-16, -3.33066907e-16, -2.22044605e-16],[ -2.22044605e-16, -3.33066907e-16, -1.94289029e-16],[ -1.11022302e-16, -3.33066907e-16, -2.49800181e-16]]]) — macrocosme, May 22 '14 at 12:07
@macrocosme That is the issue with floating point operations. You can't get values identical values after doing a bunch of operations. — Pavan Yalamanchili, May 22 '14 at 18:16
For example have a look at this: ```>>> a = math.sqrt(2) ** 2 >>> a 2.0000000000000004``` — Pavan Yalamanchili, May 22 '14 at 18:17
Working with floating point is of course always tricky! But let say one would use it to do lossless compression as for the JPEG format; it would come as a problem. — macrocosme, May 22 '14 at 23:45
@macrocosme I don't know what you are talking about. JPEG itself is a lossy compression algorithm. — Pavan Yalamanchili, May 23 '14 at 17:34

Scipy's fftpack dct and idct

1 Answers1

Linked