function B = DCT3(A)

[M, N, T] = size(A);
B = zeros(M,N,T);

for t=0:T-1
   tmpDCT2 = dct2(A(:,:,t+1));   
   for r=0:T-1
      if (r==0) 
         alpha = 1/sqrt(T);
      else 
         alpha = sqrt(2/T); 
      end 
      B(:,:,r+1) = B(:,:,r+1) + alpha * cos(pi*(2*t+1)*r/(2*T)) * tmpDCT2;
   end   
end

