function [P, G] = bld_prepost_transform(blktx, iblktx, preflt, postflt, PreOffset, PostOffset)
%Build lapped transform for the given block transform and pre/post-filters.
%The block transform can be DCT or multi-level wavelet.
%usage:
%[P, G] = bld_prepost_transform(blktx, iblktx, preflt, postflt, tile_size, PrePostLth)
%   blktx:  forward block transform,
%   iblktx: inverse block transform,
%   preflt: prefilter, size can be smaller than block size,
%   postfilter: postfilter, size can be smaller than block size, can be different from postfilter.
%           postflt may not be the inverse of the preflt.
%   PreOffset:  offset of prefilter, usually is half of the prefilter size (symmetric wrt boundary),
%   PostOffset:  offset of postfilter, usually is half of the prefilter size (symmetric wrt boundary)
%   
%   P:  rows of P are analysis filters
%   G:  this is actually the transpose of 2M x M synthesis filters.
%       The matrix reprsentation of synthesis filterbank is
%       [G1  G0   ] [yn-1]
%       [    G1 G0] [yn  ]
%       the subband inputs are filtered by [G1, G0], but the synthesis filters are defined as [G0; G1].
%   
%    PR: (G0+G1z^-1)(P0+P1z^-1) = z^-1 I.
%    ==> G0*P1 + G1 * P0 = I,   G0*P0 = G1*P1 = 0
%
%   This function can be called before, during and after FB optimization.
% Jie Liang, Simon Fraser University, Canada, jiel@sfu.ca

if nargin ~= 6
    help bld_prepost_transform;
    return;
end

tile_size = size(blktx, 1);
PreLth = size(preflt, 1);
PostLth = size(postflt, 1);

M2 = tile_size / 2; IM2 = eye(M2);  ZM2 = zeros(M2);
D0 = blkdiag(IM2, ZM2);
D1 = blkdiag(ZM2, IM2);
perm = [ZM2, IM2; IM2, ZM2];

prebdtop = tile_size / 2 - PreOffset;
Iprebdtop = eye(prebdtop);
prebdbot = tile_size - PreLth - prebdtop;
Iprebdbot = eye(prebdbot);
pre = blkdiag(Iprebdtop, preflt, Iprebdbot);

postbdtop = tile_size / 2 - PostOffset;
Ipostbdtop = eye(postbdtop);
postbdbot = tile_size - PostLth - postbdtop;
Ipostbdbot = eye(postbdbot);
post = blkdiag(Ipostbdtop, postflt, Ipostbdbot);

% fwd: (P0 + P1^z^-M) [z^-m+1; .. z^-1; 1] == > rows of [P1, P0]
P0 = blktx * D1 * perm * pre;
P1 = blktx * D0 * perm * pre;
P = [P1, P0];

% inv: [1, z-1, .. z^-m+1] * (G0 + G1^z^-M) == > col of [G0; G1]
G0 = post * perm * D0 * iblktx;
G1 = post * perm * D1 * iblktx;
G = [G0; G1]';