function demo();
% 520.643 DCT demo

% DCT transform matrix
M=8; DCT=cos(pi*[0:(M-1)]'*[1:2:(2*M-1)]/(2*M))*sqrt(2/M);DCT(1,:)=DCT(1,:)/sqrt(2);

% Verify orthogonality
DCT*DCT'
DCT'*DCT

% Plot 1D basis
bplot(DCT,DCT);

% Plot 2D basis
DCT2D(4);
DCT2D(8);

%Plot frequency responses
frplot2(DCT);

%Frequency response of the 1st basis function
plot(10*log10(abs(fft(DCT(1,:),1024))));

% Approaching KLT
Rxx=toeplitz(0.95.^[0:(M-1)]');
[KLT,D]=eig(Rxx);
KLT=flipud(KLT');
figure; bplot(KLT,DCT);
frplot2(KLT);
KLT*Rxx*KLT'
DCT*Rxx*DCT'

% DCT factorization
P=[DCT(1:2:M,:);DCT(2:2:M,:)];
U=P(1:M/2,1:M/2);
V=fliplr(P(1+M/2:M,1:M/2));
I=eye(M/2);
ButterflyM=[I fliplr(I);fliplr(I) -I];
Z=zeros(4,4);
[U Z;Z V]*ButterflyM-P

% Basis sensitivity - JND design
jnd(8,8,1);
jnd(8,8,5);
jnd(8,8,10);

% JPEG Example
x=imread('lena.bmp');
x = double(x);
figure; imshow(x/256,'truesize');
b=x(201:208,201:208)
B=DCT*b*DCT'
Bq=round(B/16)
Br=Bq*16
mydisp(b); mydisp(B); mydisp(Bq); mydisp(Br);
br=round(DCT'*Br*DCT)
e=br-b
mydisp(b); mydisp(br); mydisp(e);

b=x(17:24,201:208)
mydisp(b);
B=DCT*b*DCT'
Bq=round(B/16)
Br=Bq*16
br=round(DCT'*Br*DCT)
e=br-b
mydisp(b); mydisp(br); mydisp(e)

b=x(401:408,121:128)
mydisp(b);
B=DCT*b*DCT'
Bq=round(B/16)
Br=Bq*16
br=round(DCT'*Br*DCT)
e=br-b
mydisp(b); mydisp(br); mydisp(e)

% Blocking + Transform
fun = @dct2;
X = blkproc(x,[8 8],fun);
figure;imshow(X/max(max(X)),'truesize');
ccount(X);
X=round(X/16); 
ccount(X);
X=16*X;
xr=iglot(X,DCT); 

fun= @idct2;
xr = blkproc(X,[8 8],fun);
e=abs(x-xr);
figure; imshow(xr/max(max(xr)),'truesize');
figure; imshow(e/max(max(e)),'truesize');

% increase Q stepsize
fun = @dct2; X = blkproc(x,[8 8],fun);
ccount(X);
X=round(X/128); 
ccount(X);
X=128*X;
fun= @idct2; xr = blkproc(X,[8 8],fun);
e=abs(x-xr);
figure; imshow(xr/max(max(xr)),'truesize');
figure; imshow(e/max(max(e)),'truesize');

% Change image
x=imread('boat.tif');
x=double(x);
fun = @dct2; X = blkproc(x,[8 8],fun);
figure; imshow(X/max(max(X)),'truesize');
ccount(X);
X=round(X/16); 
ccount(X);
X=16*X;
fun= @idct2; xr = blkproc(X,[8 8],fun);
e=abs(x-xr);
figure; imshow(xr/max(max(xr)),'truesize');
figure; imshow(e/max(max(e)),'truesize');

% increase Q stepsize
fun = @dct2; X = blkproc(x,[8 8],fun);
ccount(X);
X=round(X/64); 
ccount(X);
X=64*X;
fun= @idct2; xr = blkproc(X,[8 8],fun);
e=abs(x-xr);
figure; imshow(xr/max(max(xr)),'truesize');
figure; imshow(e/max(max(e)),'truesize');

%JPEG Quantization matrix
Q=[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];

fun= @dct2; X = blkproc(x,[8 8],fun);
X=blkproc(X,[8 8],'blockQ',4);
ccount(X);
fun= @idct2; xr = blkproc(X,[8 8],fun);
e=abs(x-xr);
figure; imshow(xr/max(max(xr)),'truesize');
figure; imshow(e/max(max(e)),'truesize');

% Change image
x=imread('barb2.tif'); x=double(x);
figure; imshow(x/256,'truesize');
fun= @dct2; X = blkproc(x,[8 8],fun);
X=blkproc(X,[8 8],'blockQ',4);
ccount(X);
fun= @idct2; xr = blkproc(X,[8 8],fun);
e=abs(x-xr); 
psnr(x,xr);
figure; imshow(xr/max(max(xr)),'truesize');
figure; imshow(e/max(max(e)),'truesize');

% Wavelet
p0=[-5 0 49 0 -245 0 1225 2048 1225 0 -245 0 49 0 -5]/2048; %half-band filter
r=roots(p0); b0=[1 8 28 56 70 56 28 8 1];
q0=deconv(p0,b0);
r=[roots(q0); -ones(8,1)];
h0=poly([r(2:5)' -1 -1 -1 -1]);
f0=poly([r(1) r(6) -1 -1 -1 -1]);
h0=real(sqrt(2)*h0/sum(h0));
f0=real(sqrt(2)*f0/sum(f0));
f1=h0.*(-1).^[0:length(h0)-1];
h1=f0.*(-1).^[1:length(f0)];

X=dwt2d(x,h0,h1,3);
xmax=max(max(abs(X(1:64,1:64))))
figure; imshow(X/xmax,'truesize');
figure; imshow(X(1:64,1:64)/xmax,'notruesize');
figure; imshow(X(1:64,65:128)/max(max(abs(X(1:64,65:128)))),'notruesize');
figure; imshow(X(1:256,257:512)/max(max(abs(X(1:256,257:512)))),'notruesize');
close
figure; imshow(dwtscale(X,3)/max(max(dwtscale(X,3))),'truesize');
xr = idwt2d(X,f0,f1,3);
figure; imshow(xr/256,'truesize');
psnr(x,xr);
ccount(X);
X=32*round(X/32);
ccount(X);
close
xr=idwt2d(X,f0,f1,3);
figure; imshow(xr/256,'truesize');
psnr(x,xr);

x=imread('boat.tif'); x=double(x);
X=dwt2d(x,h0,h1,3);
figure; imshow(dwtscale(X,3)/max(max(dwtscale(X,3))),'truesize');