bindct_new.m

JPEG Image compression using IJG standards followed

function [dctcoef, bindct] = bindct_new(ver)
%Loeffler's factorization-based binDCT:
%
% Jie Liang,  May 9, 2000
%
% Note: the last butterfly is changed to a lifing,
% to keep the same number of bf in each subband.
% This adds one more shift.
%
% add cfg 8, 10 to the paper.

%Set this flag to 0 to disable the scaling of binDCT coef, used to
%check the dynamic range of the binDCt coefs.
scaling = 1;

if nargin ~= 1
   disp('Freq Response of various versions of binDCT-L:');
   disp('Usage: [dctcoef, bindct] = bindct_new(version):');
   disp('       Version: 1-->8');
   disp('       dctcoef: DCT coefficients;');
   disp('       bindct:  binDCT coefficients.');
   break;
end

M = 8;

%define the row and column indexes
ridx = [ 0 : M - 1]' ;				% column vector
cidx = ridx' ;								% row vector

%calculate the DCT coefficients
dctcoef = sqrt(2 / M) * cos(ridx * pi * (2 * cidx + 1) / (2 * M) );
dctcoef(1, : ) = dctcoef(1, : ) / sqrt(2);

close all;

liftcoef = zeros(4,3,10);

%binDCT-L1
%The dynamic range requires more bits than linDCT-C: 512 and 1024 for this case.
%Is this good or bad? more accurate approximation?
%10 shifts, 28 Adds.
%CG: 8.7716dB, Stopband: 10.307dB, MSE: 6.9201E-4.
liftcoef(:,:,1) = [1, 		1/2, 		0;		% pi/4 between X0 and X4
					   -1/2, 	1/2,		0;    % p1 and u1
                  -1/4,		1/2,		-1/4;	% p2 u2 p3
                  -1/8,		1/4,		-1/8;	% p4 u3 p5
];   

%binDCT-L2
%13 shifts, 31 Adds.
%CG: 8.8027dB, Stopband: 10.6934dB, MSE: 3.5762E-4.
liftcoef(:,:,2) = [1, 		1/2, 		0;
					   -3/8, 	1/4,		0;
                  -1/4,		1/2,		-1/4;
                  -1/8,		3/16,		-3/32;
];   

%lift (:,:,1): binDCT-New-1
%Complexity:
%16 shifts, 34 Adds.
%CG: 8.8225dB, Stopband: 9.9355dB, MSE: 4.0110E-5.
liftcoef(:,:,3) = [1, 		1/2, 		0;
					   -7/16, 	3/8,		0;
                  -1/4,	9/16,		-5/16;
                  -1/8,	3/16,		-3/32;
];   

%lift (:,:,2): binDCT-New-2
%Complexity:
%20 shifts, 38 Adds.
%CG: 8.8242dB, Stopband: 9.8926dB, MSE: 1.1179E-5.
liftcoef(:,:,4) = [1, 		1/2, 		0;
					   -13/32, 	11/32,	0;
                  -5/16,	9/16,		-5/16;
                  -3/32,	3/16,		-3/32;
];   

%lift (:,:,2): binDCT-New-2
%Complexity:
%22 shifts, 40 Adds.
%CG: 8.8257dB, Stopband: 9.9773dB, MSE: 8.1841E-6.
liftcoef(:,:,5) = [1, 		1/2, 		0;
					   -13/32, 	11/32,	0;
                  -19/64,	9/16,		-19/64;
                  -3/32,	3/16,		-3/32;
];   

%new result: 12/14/2000
liftcoef(:,:,6) = [1, 1/2, 0;
				   0, 0,  0;
                   0, 0,  0;
                   0, 0,  0;
];   

%binDCT-L1
%The dynamic range requires more bits than linDCT-C: 512 and 1024 for this case.
%Is this good or bad? more accurate approximation?
%10 shifts, 28 Adds.
%CG: 8.7716dB, Stopband: 10.307dB, MSE: 6.9201E-4.
liftcoef(:,:,7) = [1, 		1/2, 		0;		% pi/4 between X0 and X4
					   -1/2, 	1/2,		0;    % p1 and u1
                  -1/4,		1/2,		-1/4;	% p2 u2 p3
                  0,		1/4,		0;	% p4 u3 p5
];   

%7 shifts, 8.7132dB. gcd: 16, 64.
liftcoef(:,:,8) = [1, 		1/2, 		0;		% pi/4 between X0 and X4
					   -1/2, 	1/2,		0;    % p1 and u1
                  0,		1/2,		-1/4;	% p2 u2 p3
                  0,		1/4,		0;	% p4 u3 p5
];   
               
%6 shifts, 8.5623dB, GCD: 8, 32
liftcoef(:,:,9) = [1, 		1/2, 		0;		% pi/4 between X0 and X4
					   -1/2, 	1/2,		0;    % p1 and u1
                  0,		1/2,		0;	% p2 u2 p3
                  0,		1/4,		0;	% p4 u3 p5
];   


%5 shifts, 8.5464dB. gcd: 4, 16.
liftcoef(:,:,10) = [1, 		1/2, 		0;		% pi/4 between X0 and X4
					   -1/2, 	1/2,		0;    % p1 and u1
                  0,		1/2,		0;	% p2 u2 p3
                  0,		0,		0;	% p4 u3 p5
];   

%4 shifts, 8.3416dB.
liftcoef(:,:,11) = [1, 		1/2, 		0;		% pi/4 between X0 and X4
					   -1, 	1/2,		0;    % p1 and u1
                  0,		1/2,		0;	% p2 u2 p3
                  0,		0,		0;	% p4 u3 p5
];   



%Show the frequency response
freq = linspace(0, pi, 256);			% Equally spaced points between [0, 2pi].
allplot = figure;
title('DCT Frequency Response');
for i = 1 : M
   h = abs( freqz(dctcoef(i, :), 1, freq) );
%   h = h / max(h);
	figure(allplot);
   hold on;
   plot(freq, 20*log10(h));
   
   fid(i) = figure;
   plot(freq, 20*log10(abs(h)),'-.');
   v=axis;
   v(1) = 0;
   v(2) = pi;
   v(3) = -50;
   v(4) = 10;
   axis(v);
   s = sprintf('Frequency Response of DCT Basis Function No: %d', i-1);
   set(gcf, 'Name', s);
   ylabel('Magnitude Response(dB)');
   xlabel('Frequency(rad/sec)');   
end

figure(allplot);
v=axis;
v(1) = 0;
v(2) = pi;
v(3) = -50;
v(4) = 10;
axis(v);
s = sprintf('Frequency Response of All DCT Basis Functions');
set(gcf, 'Name', s);
ylabel('Magnitude Response(dB)');
xlabel('Frequency(rad/sec)');

%
% Calculate the matrix of the new BinDCT
%

%First stage: 8-row butterfly
i4 = eye(4);
j4 = fliplr(i4);
stage1 = [i4, j4;
		   j4, -i4];
      
%4-row butterfly
i2 = eye(2);
j2 = fliplr(i2);
butfly4 = [i2, j2;
		   j2, -i2];
stage2 = eye(8);
stage2(1:4,1:4) = butfly4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Upper branch rotation of pi/4: up 1; down 1/2, min x4
% for X[0] and X[4]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
predictor = eye(8);
predictor(1, 2) = liftcoef(1,1,ver);	%1

updator = eye(8);
updator(2, 1) = liftcoef(1,2,ver);		%1/2;
updator(2, 2) = -1;
rot_pi4 = updator * predictor; %without considering the scaling factors

%verify rotation pi/4
factor = eye(8);
factor(1,1) = sqrt(2)/2;
factor(2,2) = sqrt(2);
true_rot_pi4 = factor * updator * predictor; % with scaling factors

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%X[6] and X[2]:
%Upper branch rotation of 3pi/8
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
predictor = eye(8);
predictor(3, 4) = liftcoef(2,1,ver);		%-3/8; %-7/16;

updator = eye(8);
updator(4, 3) = liftcoef(2,2,ver);		%3/8; 	%11/32;

rot_pi38 = updator * predictor;

%verify rotation 3pi/8
factor = eye(8);
factor(3,3) = -sin(3*pi/8);
factor(4,4) = 1 / sin(3*pi/8);
true_rot_pi38 = factor * rot_pi38; % with scaling factors

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% rotation 3pi/16
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vlift_p1 = eye(8);
Vlift_p1(8, 5) = liftcoef(3,1,ver);

Vlift_u = eye(8);
Vlift_u(5, 8) = liftcoef(3,2,ver);

Vlift_p2 = eye(8);
Vlift_p2(8, 5) = liftcoef(3,3,ver);

low_rot_pi316 = Vlift_p2 * Vlift_u * Vlift_p1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% rotation pi/16
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vlift_p1 = eye(8);
Vlift_p1(7, 6) = liftcoef(4,1,ver);

Vlift_u = eye(8);
Vlift_u(6, 7) = liftcoef(4,2,ver);

Vlift_p2 = eye(8);
Vlift_p2(7, 6) = liftcoef(4,3,ver);

low_rot_pi16 = Vlift_p2 * Vlift_u * Vlift_p1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Two 2-row butterfly in V
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bftmp1 = [	1  0  1 0;
   			0 -1  0 1;
            1  0 -1 0;
            0  1  0 1];        
butfly2 = eye(8);
butfly2(5:8, 5:8) = bftmp1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replace the last butterfly by lifting:
%keep the same scaling of 4 in all subbands after inverse tx.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
predictor = eye(8);
predictor(8, 5) = 1;

updator = eye(8);
updator(5, 8) = 1/2;
updator(5, 5) = -1;
lastbf = updator * predictor; %without considering the scaling factors

%
%multiply and get the overall matrix,
%without consifering the final scaling parameters.
%
bindct = lastbf * butfly2 * low_rot_pi16 * low_rot_pi316 * rot_pi38 * rot_pi4 * stage2 * stage1;

%final scaling
if scaling == 1
	bindct(1,:) = bindct(1,:) * sqrt(2) / 2;
	bindct(2,:) = bindct(2,:) * (sqrt(2));
	bindct(3,:) = bindct(3,:) * (-1) *sin(3*pi/8);
	bindct(4,:) = bindct(4,:) / sin(3*pi/8);
   
   bindct(5,:) = bindct(5,:) / (sqrt(2)) * 2;   % 2 coming from the final lifting.
	bindct(6,:) = bindct(6,:) / 1;
	bindct(7,:) = bindct(7,:) / 1;
   bindct(8,:) = bindct(8,:) / (sqrt(2));
   
   bindct = bindct / 2; %inheritated from factorization.

end

%rearrange the rows
tmp = bindct;
bindct = zeros(8,8);
bindct(1,:) = tmp(1,:);
bindct(2,:) = tmp(8,:);
bindct(3,:) = tmp(4,:);
bindct(4,:) = tmp(6,:);

bindct(5,:) = tmp(2,:);
bindct(6,:) = tmp(7,:);
bindct(7,:) = tmp(3,:);
bindct(8,:) = tmp(5,:);

bplot(dctcoef, bindct);

%Show the frequency response
freq = linspace(0, pi, 256);			% Equally spaced points between [0, 2pi].
allplot = figure;
title('binDCT Frequency Response');
for i = 1 : M
   h = abs( freqz(bindct(i, :), 1, freq) );
%   h = h / max(h);
   figure(allplot);
   hold on;
   plot(freq, 20*log10(h));
   
   figure(fid(i));
   hold on;
   plot(freq, 20*log10(abs(h)));
   v=axis;
   v(1) = 0;
   v(2) = pi;
   v(3) = -50;
   v(4) = 10;
   axis(v);
   s = sprintf('Frequency Response of DCT Basis Function No: %d', i-1);
   set(gcf, 'Name', s);
   ylabel('Magnitude Response(dB)');
   xlabel('Frequency(rad/sec)');   
   legend('DCT','binDCT');
end

figure(allplot);
v=axis;
v(1) = 0;
v(2) = pi;
v(3) = -50;
v(4) = 10;
axis(v);
s = sprintf('Frequency Response of All DCT Basis Functions');
set(gcf, 'Name', s);
ylabel('Magnitude Response(dB)');
xlabel('Frequency(rad/sec)');

disp('error between DCT and binDCT:');
error = dctcoef - bindct

disp('MSE between DCT and binDCT:');
mse = sum(sum(error .* error)) / M / M;

afrplot(bindct,inv(bindct)');

Rxx=toeplitz(0.95.^[0:7]');

mse = trace( error * Rxx * error') / M

mse = trace( error * error') / M;

break;

tmp = eye(8);
tmp(6,6) = sqrt(2);
tmp(7,7) = sqrt(2);
v = tmp * butfly2 * low_rot_pi16 * low_rot_pi316;
v = v(5:8,5:8);

figure;
for r = 1 : 2
   for c = 1 : 4
      subplot(2,4,(r-1)*4+c);
	   h = abs( freqz(dctcoef((r-1)*4+c, :), 1, freq) );
      plot(freq, 20*log10(h));
      hold on;
	   h = abs( freqz(bindct((r-1)*4+c, :), 1, freq) );
      plot(freq, 20*log10(abs(h)),'-.');
      
   	v=axis;
	   v(1) = 0;
   	v(2) = pi;
	   v(3) = -50;
	   v(4) = 10;
   	axis(v);
   end
end   
subplot(2,4,5);
ylabel('Magnitude Response(dB)');
xlabel('Frequency(rad/sec)');   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% End of file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%