function [Preflt, V, codinggain, D_average, R0, R1, W, W1Nbr] = MDLTPC_design(M, N, bitrate, lostprob, cfg)
% Optimization of the prefilter and Wiener filter for the prediction compensated MD lapped transform:
% Usage: [Preflt, V, codinggain, D_average, R0, R1, W, W1Nbr] = MDLTPC_design(M, N, bitrate, lostprob, cfg)
% Input:
%   M:          Block size,
%   N:          Number of neighboring samples at each side used for Wiener filter. The size of Wiener filter is M x 2N.
%   bitrate:    Target total bit rate of two descriptions.
%   lostprob:   Probability of losing one description.
%   cfg:        0 (default): Optimize for prefilter and Wiener filter,
%               1: Fix the prefilter from the best TDLT,
%               2: Use DCT only, no prefilter.
% Output:
%   Preflt:     Optimized prefilter P,
%   V:          Optimized matrix V in the prefilter P,
%   codinggain: Coding gain of the optimized solution,
%   D_average:  Expected distortion of the optimized solution,
%   R0:         Optimized bit rate R0 for base layer,
%   R1:         Optimized bit rate R1 for enhanced layer,
%   W:          Optimized M x 2N Wiener filter,
%   W1Nbr:      Optimized M x N boundary Wiener filter for blocks that have only top/left neighboring block.
% Other output:
%   The optimized prefilter and Wiener filters are saved in MDpredictN*.dat, where the star is the value of N.
%   This file is used by the MD image coding codec.
% Examples:
%   To generate the examples in the MDLTPC paper:
%   [Preflt, V] = MDLTPC_design(8,8,1,0.2);
%   This gives V = [0.8787    0.6591    0.2426    0.1521; ... ];
%   [Preflt, V] = MDLTPC_design(8,1,1,0.2);
%   This gives V = [0.9424    0.6840    0.1942    0.1046; ... ];   
% Guoqian Sun, Jie Liang, Simon Fraser University, Canada, June 10, 2007
    
if nargin < 4
    help MDLTPC_design;
end

if nargin < 5
    cfg = 0;
end

close all;

rho = 0.95; 
 
M2 = M/2;
IM2 = eye(M2);
JM2 = fliplr(IM2);
BfM = [IM2, JM2; JM2, -IM2] / sqrt(2);
DCT = dct(eye(M));
IDCT = DCT';

if cfg == 1
    %Use best TDLT with a coding gain of 9.6151 dB
    V = [  0.9550    0.7833    0.3548    0.2391
          -0.5520    0.9008    0.6188    0.2354
           0.1124   -0.3647    1.0916    0.3904
          -0.0295    0.0082   -0.1196    1.1879];    
elseif cfg == 2
    %DCT only
    V = eye(M2);
elseif cfg == 0   
    % Optimize for prefilter and wiener filter
    % Initial value of V in prefilter
    V = eye(M2);
    ini = reshape(V,1,M2*M2);

    loops = 0;
    options = optimset('fminunc');
    myoption = optimset('MaxFunEvals',30000*length(ini), 'MaxIter',3000*length(ini),'TolX',5e-5,'TolFun',5e-5);
    options = optimset(options, myoption);      

    % Use fminunc to optimize the prefilter
    optpara = fminunc('MDLTPC_design_opt', ini, options, rho, M, N, bitrate, lostprob);

    V = reshape(optpara(1 : M2 * M2), M2, M2);
end

Preflt = BfM * blkdiag(IM2, V) * BfM;
Postflt = inv(Preflt);

% plot the frequency response
[P, G] = bld_prepost_transform(DCT, DCT', Preflt, inv(Preflt), M/2, M/2); 
bplot(P,G);
afrplot(P,G, rho);

[D_side, D_central,D_average, codinggain, R0, R1, W, W1Nbr] ...
    = MDLTPC_design_build(Preflt, rho, M, N, bitrate, lostprob);

disp(sprintf('The coding gain is %5.3f dB', codinggain));
disp(sprintf('\nR0=%f, R1=%f, R0+R1=%f', R0, R1, R0+R1));
disp(sprintf('D1 =%f, D0 =%f, D = %f',D_side, D_central, D_average));
disp(sprintf('D1D0= %g', D_side*D_central));

% Normalize Wiener filters
W =  diag(1 ./ sum(W')) * W;
W1Nbr =  diag(1 ./ sum(W1Nbr')) * W1Nbr;

%Write prefilter and Wiener filter to data file for image coding
WritePrePostFile(['MDpredictN',num2str(N),'.dat'], M, 0, 1, zeros(2*M, 2*M), Preflt, Postflt, Postflt, W, W1Nbr); 
