function Y=blkdct2(X,m,n)
% compute the 2D dct of mxn sub-blocks of X
% m must divide the # of rows of X
% n must divide the # of columns of X 
% for small blocks the dct2 is computed faster by matrix multiplication
% that by using the dct2 command
 [r,s] = size(X);
 Y = zeros(r,s);
 C1 = dct(eye(m));
 C2 = dct(eye(n))';
 for i=1:(r/m)
     for j=1:(s/n)
         e = (i-1)*m+1;
         f = (j-1)*n+1;
         Z = X(e:(e+m-1),(f:f+n-1));
         Y(e:(e+m-1),(f:f+n-1)) = C1*Z*C2;
         end;
     end;