function [normal_recon_err, bitalloc, sbnoise] = quantnoise(M, P, FwdFB, InvFB, AvgBits)
% Rqq = quantnoise(M, P, FwdFB, InvFB, AvgBits):
% Compute the quantization noise covariance matrix
% For biorthogonal FB, q noise at different subbands are different.
% Jie Liang, 06/23/2004

if nargin < 5
    help quantnoise;
    return;
end

L = size(FwdFB, 2); % columns of FB
Varx = 1; %1/3 * 128^2;
Rxx = toeplitz(0.95.^[0 : (L-1)]') * Varx;

Ryy = FwdFB * Rxx * FwdFB';
Vary = diag(Ryy);

InvFltNorm = diag(InvFB * InvFB');

WeightedNorm = Vary .* InvFltNorm;
AvgNorm = prod(WeightedNorm) ^ (1/P);

%need to multiply by the constant c in c * 2(-2b) Vary
%c = 0.75; %Gaussian signal
%Rqq = 2^(-2 * M / P * AvgBits) * AvgNorm * diag(1 ./ InvFltNorm);
%Rqq = c * 2^(-2 * M / P * AvgBits) * AvgNorm * eye(M);

%bit alloc without pos constraint
%a subband with higher WeightedNorm than AvgNorm gets more bits than average

%Note that although the total final mse = sum of (Vary_i * filternorm_i * 2^(-2*bit_i)),
%the mse of the i-th pixel is not Vary_i * filternorm_i * 2^(-2*bit_i).
%In fact, the mses of different pixels are different (reason for blocking artifact).
%the sum of (Vary_i * filternorm_i * 2^(-2*bit_i)) is obtained by using
%trace(ABC)=trace(BCA), which leads to the cond that all Vary_i * filternorm_i * 2^(-2*bit_i) should be equal.
%but the subband quant noise is Vary_i * 2^(-2*bit_i), which are not equal for biortho sys.

bitalloc = AvgBits + 1/2*log2(WeightedNorm / AvgNorm); 
sbnoise = Vary .* (2.^(-2*bitalloc));
Rqq = diag(sbnoise);

%quant error propagates to two blocks, the superimposition is the final mse for each pixel.
ReeLT = InvFB' * Rqq * InvFB;

evarLT = diag(ReeLT)';
if L == M
    normal_recon_err = evarLT;
elseif L == 2*M
    normal_recon_err = evarLT(1:M) + evarLT(M+1:2*M);
else
    disp('To be implemented.');
end