% dct2.m
% discrete cosine transform in 2d
% (painfully slow)

function F = dct2(f)

[m,n] = size(f);
F = zeros(m,n);

% transform rows
k=(2*[1:n]-1)*pi/n/2;
w=[sqrt(0.5); ones(n-1,1)]*sqrt(2/n);
for ml=1:m
   for nl=1:n
      F(ml,nl) = w(nl) * sum( f(ml,:) .* cos(k.*(nl-1)) );
   end
end

% transform columns
f=F;
k=(2*[1:m]'-1)*pi/m/2;
w=[sqrt(0.5); ones(m-1,1)]*sqrt(2/m);
for nl=1:n
   for ml=1:m
      F(ml,nl) = w(ml) * sum( f(:,nl) .* cos(k.*(ml-1)) );
   end
end