function [gdb,Ryy, opt_recerr, bitalloc, sbnoise] = bgtc(H, G, bitalloc, bitbudget, rho);
% [gdb,Ryy, opt_recerr, bitalloc, sbnoise]=bgtc(H,G);
% Function to calculate the coding gain of a biorthogonal PRFB, does not apply to non-PR FB.
% Rows of H hold the analysis filter coefficients
% Rows of G hold the synthesis filter coefficients (G is the transpose of the inverse transform)
% bitalloc: 0 no positive constraint,
%           1 use positive constraint
%   bitbudget:  average bit used in bitalloc = 1.
% opt_recerr_orgorder: 
%       Optimal output noise in the reconstructed signal caused by each subband.
%-----------------------------------------------------------------------

if nargin == 2
    bitalloc = 0;
    % The bits will be canceled in the classcal coding gain def
    % so any number here is ok.
    bitbudget = 8;  
    rho = 0.95;
end

[m,n]=size(H);
Rxx=toeplitz(rho.^[0 : (n-1)]');
Ryy=H*Rxx*H';
s=diag(Ryy)';
% g=mean(s)/(prod(s)^(1/m));
% g=1/((prod(s.*sqrt(diag(G*G'))))^(1/m));

Gnorm = diag(G*G')';
if find(Gnorm < 1E-10)
    disp('Zero column exists in inverse transform.');
    return;
else
    if bitalloc == 0
        %allow neg bit allocation, use classic bit allocation and coding gain def
        %opt_recerr is the distribution of recon error at each pixel
        [opt_recerr, bitalloc, sbnoise] = quantnoise(m, m, H, G, bitbudget);
    else
        %pos bits
        [opt_recerr, bitalloc, sbnoise] = posbitalloc(H, G, bitbudget);           
    end
    
    %if neg bits allowed, the bit term will be canceled, this reduces to the classc coding gain def
    gdb = m * 2^(-2 * bitbudget) / sum(opt_recerr);     
    
    gdb=10*log10(gdb);
end