bindct_c.m

binDCT fron D.Tran used in Image compression

function [DCT, binDCT, binDCT4, scoef]=binDCT_C(ver, scale, toplot, liftbf)
%binDCT Type C (Chen-Wang's algorithm):
%usage: 
% [DCT, binDCT, binDCT4, s]=binDCT_C(ver, scale, plot, liftbf)
%     ver:      1 to 9, as in the IEEE SP paper.
%     scale:    1: retun scaled binDCT. 0: not scale.
%     toplot:   1: plot frequency response.
%     liftbf:   0: use butterflies, 1: replace all butterflies by lifting.
% return:
%     DCT:     DCT coefficient;
%     binDCT:  binDCT coefficent;
%     binDCT4: 4-point binDCT.
%     scoef:   Scaling vector in natural order.
%
% Reference:
% J. Liang, T. D. Tran, Fast Multiplierless Approximations of the DCT with the Lifting
% Scheme, IEEE Trans. Signal Processing, Vol. 49, No. 12, pp. 3032-3044, Dec. 2001. 
%
% Trac D. Tran and Jie Liang
% ECE Department, The Johns Hopkins University
% 3400 North Charles Street, 105 Barton Hall,
% Baltimore, MD 21218
% E-mail: [e-mail], [e-mail]
% Dec. 2000
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Copyright (c) 2000 Trac D Tran and Jie Liang
% This program is Copyright (c) by Trac D Tran and Jie Liang.
% It may not be redistributed without the consent of the copyright
% holders. In no circumstances may the copyright notice be removed.
% The program may not be sold for profit nor may they be incorporated
% in commercial programs without the written permission of the copyright
% holders. This program is provided as is, without any express or
% implied warranty, without even the warranty of fitness for a
% particular purpose.
%-----------------------------------------------------------------------

if nargin ~= 4
   help binDCT_C;
   break;
end   

DCT = dct(eye(8));

%		   p1, 	u1, 	p2, 	u2, 	p3, 	u3,	    p4, 	u4, 	p5
para = [	...
      	13/32	11/32	11/16	15/32	3/16	3/16	13/32	11/16	13/32; ...
         7/16	3/8	    5/8	    7/16	3/16	3/16	13/32	11/16	13/32; ...
         13/32	11/32	11/16	15/32	3/16	3/16	7/16	11/16	3/8; ...
         7/16	3/8	    5/8	    7/16	3/16	3/16	7/16	11/16	3/8; ...
      	3/8	    3/8	    7/8	    1/2	    3/16	3/16	7/16	11/16	3/8; ...
      	1/2	    3/8	    7/8	    1/2	    3/16	1/4	    7/16	3/4	    3/8; ...
         1/2, 	1/2, 	1, 	    1/2,	1/4,	1/4, 	1/2,	3/4, 	1/2; ...
         1  	1/2	    -1		-1/2	    0		0	    0	    1/2		1/2; ... 
         0		0		0		0		0		0		0		0		0; ...
      ];
  

%norm for butterflies
if liftbf == 0
   norm = 1;
else
   norm = 2;
end

I4 = eye(4);	J4 = fliplr(I4);
bf1 = [I4/norm, J4/norm; J4, -I4];

%matrix U
I2 = eye(2);	J2 = fliplr(I2);
bf2 = [I2/norm, J2/norm; J2, -I2];
stage3 = blkdiag([1 1; 1/2 -1/2], [1 0; -para(ver,2), 1]*[-1 para(ver,1); 0 1]);
U = stage3 * bf2;

%matrix V
V1 = eye(4);
V1(2:3,2:3) = [-1, para(ver,9); 0 1] * [1 0; para(ver,8), 1] * [1 -para(ver,7); 0 1];

V20 = [1/norm 1/norm; 1 -1];
V21 = [-1 1; 1/norm 1/norm];
V2 = blkdiag(V20, V21);

V3 = eye(4);
V3([2,3],[2,3]) = [1 0; -para(ver,4), 1] * [1 para(ver,3); 0 1];
V3([1,4],[1,4]) = [1 0; -para(ver,6), 1] * [-1 para(ver,5); 0 1];

V = V3 * V2 * V1;

binDCT = blkdiag(U, V) * bf1;

%rearrangement
tmp = binDCT;
bindct = zeros(8,8);
binDCT(1,:) = tmp(1,:);
binDCT(5,:) = tmp(2,:);
binDCT(7,:) = tmp(3,:);
binDCT(3,:) = tmp(4,:); 
binDCT(8,:) = tmp(5,:);
binDCT(6,:) = tmp(6,:);
binDCT(4,:) = tmp(7,:);
binDCT(2,:) = tmp(8,:);

if scale == 1
   %in natural order: for X[0] to X[7]:
   scoef = [sin(pi/4)/2* (norm^2), 	1/(2*sin(7*pi/16))*norm, 	1/2/sin(3*pi/8)*norm,	1/(2*cos(3*pi/16)), ...
        sin(pi/4)*(norm^2), 		cos(3*pi/16)/2, 				sin(3*pi/8)/2*norm, 	sin(7*pi/16)/2*norm];
   binDCT = diag(scoef) * binDCT;
end

if toplot == 1
    bplot(DCT, binDCT);
    set(gcf, 'name','Basis Functions of 8-point DCT (left) and binDCT (right)');        
	afrplot(DCT, DCT');
    set(gcf, 'name','Frequency Response of 8-point DCT');    
	afrplot(binDCT, inv(binDCT)');
    set(gcf, 'name','Frequency Response of 8-point binDCT');
end
gain = bgtc(binDCT, inv(binDCT)');
s=sprintf('8-point binDCT coding gain: %fdB.',gain);
disp(s);

%%%%%%%%%%%%%%%%%%
% 4-point binDCT: from U
%%%%%%%%%%%%%%%%%%
DCT4 = dct(eye(4));
binDCT4 = zeros(4,4);
binDCT4 (1,:) = U(1,:);
binDCT4 (3,:) = U(2,:);
binDCT4 (4,:) = U(3,:);
binDCT4 (2,:) = U(4,:);

if scale == 1
   s = [1/2 * norm, 1/(sqrt(2)*sin(3*pi/8)), 1* norm, sin(3*pi/8)/sqrt(2)];
   binDCT4 = diag(s) * binDCT4;
end   

if toplot == 1
    bplot(DCT4, binDCT4);
    set(gcf, 'name','Basis Functions of 4-point DCT (left) and binDCT (right)');    
	afrplot(DCT4, DCT4');
    set(gcf, 'name','Frequency Response of 4-point DCT');    
	afrplot(binDCT4, inv(binDCT4)');
    set(gcf, 'name','Frequency Response of 4-point binDCT');    
end
gain = bgtc(binDCT4, inv(binDCT4)');
s=sprintf('4-point binDCT coding gain: %fdB.',gain);
disp(s);

% Calculate the MSE
if scale == 1
    error = binDCT - DCT;
    Rxx=toeplitz(0.95.^[0:7]');
    mse = trace( error * Rxx * error') / 8;
    s=sprintf('MSE between DCT and binDCT: %0.4g\n', mse);
    disp(s);
end