clear all;
Del = [2 4 8 10; 4 4 8 10; 8 8 8 10; 10 10 10 10];
T = 0.5*[1 1 1 1 ; 1 1 -1 -1; 1 -1 -1 1; 1 -1 1 -1];
f = [ 50 50 70 70 ; 50 50 80 80 ; 40 40 90 90 ; 40 40 90 90 ];
%transform coefficients
fbv = T*(f-50);
fb = fbv * T';
fbd = fb./Del;
%round down
eps = .0000001;
%output of quantizer encoders
fbe = round(fb./Del-eps*sign(fb));
%output of quantizer decoders
fbd = fbe.*Del;
%reconstructed image
fhat = T'*(fbe.*Del)*T+50;
eb = fb-fbd;
e = f-fhat;
%Mean-square in the transform domain
MSEt = sum(sum(eb.*eb))/16;
%Mean-square error in the image domain
MSE = sum(sum(e.*e))/16;
PSNR = 10*log10(100^2/MSE);