function G = CodingGain_num(h0_coeff, h1_coeff, r, level, model, M)
% compute the coding gain G_{SBC} for octave-band quincunx filter banks
% input:  h0_coeff, h1_coeff -- the pair of analysis filters coefficients in matrix form
%                          r -- AR model parameter (usually taken as 0.95)
%                      level -- levels of decomposition
%                      model -- 1; seperable model; 0: isotropic model
%                          M -- sampling matrix
% output:                  G -- a vector of the coding gains for increasing levels of decompostion
% e.g. G = CodingGain_num(h0_coeff, h1_coeff, 0.95, 6, 0, [1 1; 1 -1])
% Copyright (c) 2006 Yi Chen

t = cputime;
[h0_coeff, td0] = reduce_size(h0_coeff, [0, 0], 1e-10);
[h1_coeff, td1] = reduce_size(h1_coeff, [0, 0], 1e-10);

% compute the synthesis filter coefficients from the analysis pair
g0_coeff = h1_coeff;
for i = 1:length(g0_coeff(:,1))
    j = 2-mod(i,2):2:length(g0_coeff(1,:));
    g0_coeff(i,j) = -g0_coeff(i,j);
end
g1_coeff = h0_coeff;
for i = 1:length(g1_coeff(:,1))
    j = 2-mod(i,2):2:length(g1_coeff(1,:));
    g1_coeff(i,j) = -g1_coeff(i,j);
end

A1 = zeros(1,level);   % A1(k) for analysis highpass channel in level k
B1 = zeros(1,level);   % B1(k) for synthesis highpass channel in level k
A0 = 0; % analysis lowpass
B0 = 0; % synthesis lowpass
G = zeros(1,level);   % coding gain of level k
[h1_coeff, td] = reduce_size(h1_coeff, [0, 0], 1e-10); 
[g1_coeff, td] = reduce_size(g1_coeff, [0, 0], 1e-10);
[h0_coeff, td] = reduce_size(h0_coeff, [0, 0], 1e-10);
[g0_coeff, td] = reduce_size(g0_coeff, [0, 0], 1e-10);
h1 = h1_coeff; 
g1 = g1_coeff;
h0 = h0_coeff;
g0 = g0_coeff;

A0 = calc_A(h0, r, model);
A1(1) = calc_A(h1, r, model);
B0 = sum(sum(g0.^2));
B1(1) = sum(sum(g1.^2));

Gt = (A1(1)*B1(1))^(1/2);    
G(1) = 10*log10(1/(Gt*(A0*B0)^(1/2)));

for i = 2:level
    [h0t, td] = upsampling_2d(h0, [0,0], M);
    h0 = conv2(h0_coeff,h0t);
    [h0, td] = reduce_size(h0,  [0,0], 1e-5);
    A0 = calc_A(h0, r, model);
    
    [h1t, td] = upsampling_2d(h1, [0,0], M);
    h1 = conv2(h0_coeff,h1t);
    [h1, td] = reduce_size(h1, [0, 0], 1e-5);
    A1(i) = calc_A(h1, r, model);
     
    [g0t, td] = upsampling_2d(g0, [0,0], M);
    g0 = conv2(g0_coeff,g0t);
    B0 = sum(sum(g0.^2));
    
    [g1t, td] = upsampling_2d(g1, [0,0], M);
    g1 = conv2(g0_coeff,g1t);
    B1(i) = sum(sum(g1.^2));
    
    Gt = Gt * (A1(i)*B1(i))^(1/2^i);
    G(i) = 10*log10(1/(Gt*(A0*B0)^(1/2^i)));
end    

t = cputime-t;

function A = calc_A(h, r, model)    
% local function to calculate A_k with filter coefficients h and AR model
% parameter r

m = size(h, 1);
n = size(h, 2);

M = zeros(m, m, n);
if model == 0
    for i = 0:n-1
        a = r.^(sqrt((0:m-1).^2 + i^2));
        M(:, :, i+1) = toeplitz_c(a);
    end
else
    for i = 0:n-1
        a = r.^((0:m-1) + i);
        M(:, :, i+1) = toeplitz_c(a);
    end
end

hv = reshape(h, [numel(h), 1]);
A = 0;
for i = 1:n
    Mt = reshape(M(:, :, abs(i-(1:n))+1), [m, m*n]);
    A = A + h(:, i)'*(Mt*hv);
end

clear M*