function [D_side, D_central, D_average, codinggain, R0, R1, W, W1Nbr] =...
    MDLTPC_design_build(Preflt, rho, M, N, bitrate, lostprob)
% Routine called by MDLTPC_design.m
% Guoqian Sun, Jie Liang, Simon Fraser University, Canada, June 10, 2007

% Get DCT, prefilter, postfilter etc
R = bitrate;

C = dct(eye(M));
C33 = blkdiag(C, C, C);
C22 = blkdiag(C, C);

Pre0 = Preflt(1 : M/2, :);
Pre1 = Preflt(M/2+1 : M, :);
Pre34 = blkdiag(Pre1, Preflt, Preflt, Pre0);

Postflt = inv(Preflt);
Post0 = Postflt(:,1:M/2);
Post1 = Postflt(:,M/2+1:M);
Post21 = blkdiag(Post1, Post0);

% Build the lapped transform, both P and G are M x 2M
[P, G] = bld_prepost_transform(C, C', Preflt, inv(Preflt), M/2, M/2); %forward and backward transform

% Build R_{S_3,S_3}
Rx4x4 =  toeplitz(rho .^ [0 : (4 * M-1)]'); 
Rs123s123 = Pre34 * Rx4x4 * Pre34';

% Get submatrices for Wiener filter
Rs13s13 = [Rs123s123(1:M,1:M),       Rs123s123(1:M,2*M+1:3*M);
           Rs123s123(2*M+1:3*M,1:M), Rs123s123(2*M+1:3*M,2*M+1:3*M)];
Rs2s13 = [Rs123s123(M+1:2*M,1:M),    Rs123s123(M+1:2*M,2*M+1:3*M)];

% R_{\hat{S}_3, \hat{S}_3}: Ignoring quantization error C33' * blkdiag(Kqq, Kqq, Kqq) * C33;
Rs123s123_hat = Rs123s123;  

Rs13s13_hat =  [Rs123s123_hat(1:M, 1:M),  Rs123s123_hat(1:M,      2*M+1:3*M);
          Rs123s123_hat(2*M+1:3*M, 1:M), Rs123s123_hat(2*M+1:3*M, 2*M+1:3*M)];
       
% Generate the M x 2N Wiener filters
Rss = Rs123s123(1:M,1:M);
Rs13Ls13L_hat = Rs13s13_hat(M-N+1 : M+N, M-N+1 : M+N);  % 2N x 2N
W = Rs2s13(:, M-N+1 : M+N) * inv(Rs13Ls13L_hat);       % M x 2N Wiener filter
W = [zeros(M,M-N), W, zeros(M,M-N)];    % Zero padding to M x 2M, to simplify the implementation.

% Boundary M x N Wiener filter, for blocks that have only top/left neighboring block.
% The filter for blocks with bottom/right neighbor is: fliplr(flipud(W1Nbr));  
Rs1N = Rs13s13_hat(M-N+1 : M, M-N+1 : M);
W1Nbr = Rs2s13(:, M-N+1 : M) * inv(Rs1N);
W1Nbr = [zeros(M,M-N), W1Nbr];

Rs2s2 = Rs123s123(M+1:2*M, M+1:2*M);
Rdd = Rs2s2 - W * Rs2s13';   % Min value of Rdd when Wiener filter is used.
Ruu = C * Rdd * C';

Ryy = C * Rss * C';
sigma_y = diag(Ryy);
sigma_d = diag(Rdd);
sigma_u = diag(Ruu);

Gnorm = diag(G*G')';    % G: M x 2M.

%Geometric mean: sigma^2_{B,M} and sigma^2_{E,M} for base layer and enhanced layer, respectively
sigma_0 = prod(sigma_y .* Gnorm') ^ (1/M);
sigma_1 = prod(sigma_u .* Gnorm') ^ (1/M);

% Bit allocation
if (lostprob == 0)
    delta = inf;
else
    delta = (1/4) * log2(sigma_0 / (lostprob * sigma_1));
end
R0 = R/2 + delta;
R1 = R/2 - delta;
R0 = min(R, R0);
R1 = max(0, R1);

% D_{B,M} and D_{E,M}
epsilon = 1;   
D_MB = epsilon * sigma_0 * 2^(-2 * R0);
D_ME = epsilon * sigma_1 * 2^(-2 * R1);

% Side distortiona nd central distortion
D_side = (D_MB + D_ME) / 2;
D_central = D_MB;

%Expected distortion
D_average = (1-lostprob)^2 * D_central + 2 * lostprob * (1-lostprob) * D_side;

codinggain = bgtc(P, G, 0, bitrate, rho);

% End