function non_errorfree_shannon_limit()
% The key equation is:
% R*(1-H(Pe)) = C(R*Eb/N0)
clear; clc; close all;
R = .5;
Pe_array = [2*1e-3, 10.^(-[3:.25:6])];
EbN0_array = Pe_array*0;
for n = 1:length(Pe_array)
    Pe = Pe_array(n);
    T = [];
    u = -0.01:5e-4:0.2;
    for EbN0 = u
        t = R*(1-entropy(Pe)) - C(R, EbN0);
        T = [T, abs(t)];
    end
    EbN0_array(n) = u(T == min(T));
end
semilogy(EbN0_array, Pe_array,'+-');
grid on; xlabel('E_b/N_0'); ylabel('P_e');


function capacity = C(R, EbN0)
ebn0 = 10^(EbN0/10);
sigma2 = 1/2/R/ebn0;
integral_delta = sqrt(sigma2)/100;
y = -5*sqrt(sigma2)-1 : integral_delta : 5*sqrt(sigma2)+1;
f_0 = 1 / sqrt(2*pi*sigma2) * exp(-(y-1).^2./(2*sigma2)) + eps;
f_1 = 1 / sqrt(2*pi*sigma2) * exp(-(y+1).^2./(2*sigma2)) + eps;
c0 =  sum(f_0 .* log2(2*f_0./(f_0+f_1))) * integral_delta;
c1 =  sum(f_1 .* log2(2*f_1./(f_0+f_1))) * integral_delta;
capacity = (c0+c1)/2;

function y = entropy(x)
y = x*log2(x+eps) + (1-x)*log2(1-x+eps);
y = -y;