function [Dist, Rate,B_idct] = blk_calc(B, step)


  %% PER BIT VALUES FOR DISTORTION & RATE

   M = size(B, 2);

   B_dct = dct2(B);
   B_q = quant(B_dct, step); 
   B_idct = idct2(B_q); 

   % Find the average distortion
   Dist = sum(sum((B_q - B_dct).^2))/(M*M); 

   % Find the B_q bitrate
   range = max(max(B_q)) - min(min(B_q));
   pmf = zeros(range, 1);
   temp_pmf = hist(B_q, range);
   for n=1:size(temp_pmf,2)
      pmf = temp_pmf(:,n) + pmf;
   end
   pmfsize = size(pmf);
   pmf = pmf/sum(pmf);

   % bits/pixel*M*M = bits/B_q
   Rate = sum( -pmf.*log2(pmf + (pmf ==0)) );