% Usage: [x2,P] = jpegdemo(x,r)
% Simple JPEG compression demo. Input argument x is an M by N
% array of numbers, either uint8 or double (converted to double
% below). Argument "r" is for scaling the quantization matrix as
% in book. Uses quantization matrix from book. Output is M by N
% array "x2", simulated compressed/decompressed image as an array
% of doubles, and also "P", fraction of non-zero block DCT
% coefficients (as a crude measure of compression). Calls the routine
% blkidct2().
% If you don't have the signal processing toolbox substitute kdct()
% for dct().
function [x2,P] = jpegdemo(x,r)
T = r*[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 194 113 92;
49 64 78 87 103 121 120 101; 72 92 95 98 121 100 103 99];
m=8; n=8; %Really, the blocksize m,n can be anything.
% Extend image x so dimensions divisible by m, n.
[M0,N0] = size(x); M=M0; N=N0;
blocksM = ceil(M/m); blocksN = ceil(N/n);
M = m*blocksM; N = n*blocksN;
for k=(M0+1):M %Extend each column
x(k,:) = x(M0,:);
end
for k=(N0+1):N %Extend each row
x(:,k) = x(:,N0);
end
Y = zeros(M,N); %Storage for block DCT
%Arrays for 8x8 DCT; straight matrix multiply is faster than dct2.
C1 = dct(eye(m));
C2 = dct(eye(n))';
%Block DCT code, quantization
for i=1:blocksM
for j=1:blocksN
%Pick off each 8 x 8 block, put in "Z"
e = (i-1)*m+1;
f = (j-1)*n+1;
Z = double(x(e:(e+m-1),(f:f+n-1)));
%Compute DCT of block
Q = C1*Z*C2;
%Quantize (quant. matrix already scaled), install in Y.
Q2 = T.*round(Q./T);
Y(e:(e+m-1),(f:f+n-1)) = Q2;
end;
end;
P = sum(sum(sign(abs(Y))))/M/N;
%Inverse block DCT
Y = double(blkidct2(Y,m,n));
%Pick off original image dimensions, return as a double
x2 = Y(1:M0,1:N0);