How to calculate FDCT

Question

I'm new to Mathemtica and I'm trying to calculate Discrete Cosine Transformation FDCT. I found the FourierDCT built-in function, but not DCT, so I need to implement it. I have tried couple of ideas but no luck yet. Any ideas that can get me started will be really appreciated.

Edit 1 ******************************

This is the result of FourierDCT on the same Matrix from the Wikipedia link above, and the values are different. I need Forward Discrete Cosine Trasnformation, not the FourierDCT:

From Wikipedia: "A DCT is similar to a Fourier transform in the sense that it produces a kind of spatial frequency spectrum."

Edit 2 ******************************

This is the matrix I'm using as input in Mathemtica:

b = {{52, 55, 61, 66, 70, 61, 64, 73}, {63, 59, 55, 90, 109, 85, 69, 72}, {62, 59, 68, 113, 144, 104, 66, 73}, {63, 58, 71, 122, 154, 106, 70, 69}, {67, 61, 68, 104, 126, 88, 68, 70}, {79, 65, 60, 70, 77, 68, 58, 75}, {85, 71, 64, 59, 55, 61, 65, 83}, {87, 79, 69, 68, 65, 76, 78, 94} }

Answer 1

The cosines of frequencies different to zero, have half the squared area than unity, so 1/Sqrt@2 times the norm. So, the non-first rows and columns need to be scaled up in a factor of Sqrt@2 (smaller dual basis, higher basis, smaller coordinates, need to increase them to normalize). At least this seems to fit the sample data you gave

normalizedDCT[b_] := Module[{res = 2` FourierDCT[b - 128]},
 res[[1]] = res[[1]]/Sqrt@2`;
 res[[All, 1]] = res[[All, 1]]/Sqrt@2`;
 res
 ]

normalizedDCT[b]//MatrixForm

Shorter, but a tiny bit slower

normalizedDCT[b_] := 
     With[{d = DiagonalMatrix@SparseArray[{1 -> 1`}, Length@b, Sqrt[2`]]},
         d . FourierDCT[b - 128].d]

Answer 2

I do not know how to incorporate that scaling factor shown in the JPEG wiki page into the DCT formula used by Mathematica. But here is a quick hack implementation to implement that DCT formula used in the jpeg wiki page you linked to. This gives the same matrix now as shown there. Feel free to make it more 'functional' (i.e. not using For loops).

g[u_, v_, b_] := Module[{val = 0, alpha1, alpha2},       
   alpha1 = If[u == 0, Sqrt[1/8], Sqrt[2/8]];
   alpha2 = If[v == 0, Sqrt[1/8], Sqrt[2/8]];

   For[x = 0, x <= 7, x++,
    For[y = 0, y <= 7, y++,
     val +=b[[x + 1, y + 1]]*alpha1*alpha2*Cos[Pi/8 (x + 1/2)*u]*Cos[Pi/8 (y + 1/2)*v]
     ]
    ];

   val
   ];

b = {{52, 55, 61, 66, 70, 61, 64, 73}, {63, 59, 55, 90, 109, 85, 69, 
    72}, {62, 59, 68, 113, 144, 104, 66, 73}, {63, 58, 71, 122, 154, 
    106, 70, 69}, {67, 61, 68, 104, 126, 88, 68, 70}, {79, 65, 60, 70,
     77, 68, 58, 75}, {85, 71, 64, 59, 55, 61, 65, 83}, {87, 79, 69, 
    68, 65, 76, 78, 94}};

b = b - 128;

Table[g[u, v, b], {u, 0, 7}, {v, 0, 7}] // N // MatrixForm

The result is

which is the same as wiki page.