bindct_l.m

binDCT fron D.Tran used in Image compression

function [DCT, binDCT, scoef]=binDCT_L(ver, scale, toplot, liftbf)
%binDCT Type L (Loeffler's algorithm):
%usage: 
% [DCT, binDCT, s]=binDCT_L(ver, scale, plot)
%     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;
%     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_L;
    return;   %break; 
end   

DCT = dct(eye(8));

%		   p1, 	u1, 	p2, 	u2, 	p3, 	p4, 	u3, 	p5
para = [	...
      	13/32	11/32	19/64,  9/16,   19/64,  3/32,   3/16,   3/32; ...
      	13/32	11/32	5/16,   9/16,   5/16,   3/32,   3/16,   3/32; ...
      	7/16	3/8 	1/4,    9/16,   5/16,   1/8,    3/16,   3/32; ...
      	3/8 	1/4 	1/4,    1/2,    1/4,    1/8,    3/16,   3/32; ...
      	1/2 	1/2 	1/4,    1/2,    1/4,    1/8,    1/4,    1/8; ...        
      	1/2 	1/2 	0,      1/2,    1/4     0       1/4     0; ...
      	1/2 	1/2 	0,      1/2,    0       0       0       0; ...
      	1    	1/2 	0,      1/2,    0       0       0       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([1,4],[1,4]) = [1, 0; -para(ver,5), 1] * [1 para(ver,4); 0, 1] * [1, 0; -para(ver,3), 1];

V2 = eye(4);
V2([2,3],[2,3]) =  [1, 0; -para(ver,8), 1] * [1 para(ver,7); 0, 1] * [1, 0; -para(ver,6), 1];

V3 = eye(4);
V3([1,3], [1,3]) = [1/norm 1/norm; 1 -1];
V3([2,4], [2,4]) = [-1 1; 1/norm 1/norm];

V4 = eye(4);
V4([1,4],[1,4]) = [-1, 1/2; 0, 1] * [1, 0; 1, 1];

V = V4 * 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(4,:) = tmp(6,:);
binDCT(6,:) = 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/sqrt(8)*norm, 	1/sin(3*pi/8)/2 * norm,	1/2, ...
        sin(pi/4)*(norm^2), 		1/2, 				sin(3*pi/8)/2*norm, 	1/sqrt(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);

% 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