%Study of DCT quantization using JPEG quantization matrix
clear all
b = im2double(imread('barb.tif'));
npix = prod(size(b));
T = dctmtx(8);
BB = blkproc(b, [8 8], 'P1*x*P2', T, T'); %compute block DCT
k=0;
for qf = .1:.1:3
   k=k+1;
   Del = qf*[16 11 10 16 24 40 51 61;
   12 12 14 19 26 58 60 55;
   14 13 16 24 40 57 69 56;
   14 17 22 29 51 87 80 62;
   18 22 37 56 68 109 103 77;
   24 35 55 64 81 104 113 92;
   49 64 78 87 103 121 120 101;
   72 92 95 98 112 100 103 99]'/255;
   
   BBhat = blkproc(BB, [8 8], 'round(x./P1).*P1',Del);%quantize with JPEG matrix
   bhat = blkproc(BBhat, [8 8], 'P1*x*P2', T', T); %inverse block DCT
   %Compute PSNR
   MSE = sum(sum((b-bhat).^2))/npix
   PSNR(k) = 10*log10(1/MSE);
end
qflist = .1:.1:3;
plot(qflist,PSNR)
xlabel('qf')
ylabel('PSNR')
title('PSNR as a function of parameter qf for JPEG quantization')