function [D0, D1, D, PD0, PD1, PD, R] = MDLTPC_DPCM(rho, pe)
% Compare the performance of Jayant's method and prediction compensated MDC in the DPCM case.
% This is Fig. 3 of the MDLTPC paper.
% The product PD0(31) * PD1(31) gives the D0D1 values of DPCM in Table I: 1.88, 1.95, 2.04, 2.23.
% Other numbers in Table I are given by MDLTPC_design.m.
% Usage:
%   Fig3_MDLTPC_DPCM(rho, pe)
%       rho: correlation coefficient
%       pe: probabilty of losing one description
% Example:
%   MDLTPC_DPCM(0.95, 0.1);
%
% Note: R(i) is bit/sample of ALL descriptions.
% Jie Liang, Simon Fraser University, Canada, jiel@sfu.ca, July 2008

if nargin < 2
    help MDLTPC_DPCM;
    return;
end

p0 = ((1-pe)^2);
p1 = 2*pe*(1-pe);

epsilon = 1;

Rstep = 0.1;
R = 1 : Rstep : 6;
lth = length(R);

D0 = zeros(lth, 1);
D1 = D0;
D = D0;

% Jayant's DPCM, no coding of pred error
for i = 1 : lth;
    D1(i) = 1/2 * epsilon * ((1-rho^4)*2^(-2*R(i)) + (1-rho^2)/(1+rho^2));
    D0(i) =  epsilon * (1-rho^4)*2^(-2 * R(i));
    D(i) = p0 * D0(i) + p1 * D1(i);
end

%Prediction-compensated MDC: need opt bit alloc
PD0 = zeros(lth, 1);
PD1 = PD0;
PD = PD0;
for i = 1 : lth
    R0(i) = 1/2*R(i) + 1/4*log2((1+rho^2)^2/pe);
    R1(i) = 1/2*R(i) - 1/4*log2((1+rho^2)^2/pe);
    
    R0(i) = min([R(i), R0(i)]);
    R1(i) = max([0, R1(i)]);
    
    PD1(i) = 1/2 * epsilon * ((1-rho^4)*2^(-2*R0(i)) + (1-rho^2)/(1+rho^2)*2^(-2*R1(i)));
    PD0(i) =  epsilon * (1-rho^4)*2^(-2 * R0(i));
    PD(i) = p0 * PD0(i) + p1 * PD1(i);
end


figure;
hold on; plot(R(1:10:lth), 10*log10(1./D0(1:10:lth)), 'rd', 'LineWidth', 2);
hold on; plot(R(1:10:lth), 10*log10(1./PD0(1:10:lth)),'r^', 'LineWidth', 2);
hold on; plot(R(1:10:lth), 10*log10(1./D(1:10:lth)), 'k+', 'LineWidth', 2);
hold on; plot(R(1:10:lth), 10*log10(1./PD(1:10:lth)), 'ks', 'LineWidth', 2);
hold on; plot(R(1:10:lth), 10*log10(1./D1(1:10:lth)), 'bv', 'LineWidth', 2);
hold on; plot(R(1:10:lth), 10*log10(1./PD1(1:10:lth)),'bo', 'LineWidth', 2);

legend('1/D_0: Jayant', '1/D_0: Proposed', '1/D: Jayant', '1/D: Proposed', '1/D_1: Jayant', '1/D_1: Proposed');
xlabel('Total Rate (bits/sample)'); ylabel('SNR (dB)');
grid on;

hold on; plot(R(1:lth), 10*log10(1./D0(1:lth)), 'r:', 'LineWidth', 2);
hold on; plot(R(1:lth), 10*log10(1./PD0(1:lth)),'r-', 'LineWidth', 2);
hold on; plot(R(1:lth), 10*log10(1./D(1:lth)), 'k:', 'LineWidth', 2);
hold on; plot(R(1:lth), 10*log10(1./PD(1:lth)), 'k-', 'LineWidth', 2);
hold on; plot(R(1:lth), 10*log10(1./D1(1:lth)), 'b:', 'LineWidth', 2);
hold on; plot(R(1:lth), 10*log10(1./PD1(1:lth)),'b-', 'LineWidth', 2);

return;
