function [lower] = BINCAPLO(SNR, capacity)
% BINCAPLO - Lower bound on binary input AWGN channel capacity
%   C = BINCAPLO(SNR) calculates a lower bound C on the capacity
%   of the binary input AWGN channel for a given SNR (in dB).

% Note: Due to problems with numerical precision, the original bounds gave 
% inaccurate results at high (>12 dB) SNR.  As a temporary fix, for high SNR values 
% we use the bounds with fewer numbers of Taylor series terms.  Work is currently 
% underway to reimplement the bounds in a more numerically stable manner.

%  (C) Mark Shane (shane@ee.ucla.edu) and Richard Wesel (wesel@ee.ucla.edu)
%  University of California, Los Angeles

if nargin == 1, capacity = 0; end

lowSNR = find(SNR < 12);
hiSNR  = find(SNR >= 12);

lowsnr = 10.^(SNR(lowSNR)/10);
hisnr = 10.^(SNR(hiSNR)/10);

snr = lowsnr;
a = 1.1906./snr;
d_taylow4 = - log(2) - snr.*(1-Q(-sqrt(snr)) + Q(sqrt(snr))) ...
	  + sqrt(2*snr/pi).*exp(-snr/2) + (Q((a+1).*sqrt(snr)) + ...
     exp(4*snr).*Q((a+3).*sqrt(snr))) ...  % original bound
     ... 
     + log(2)*(1-Q((a-1).*sqrt(snr)) - Q((a+1).*sqrt(snr))) ... % ln(2)
     ...
     + sqrt(snr/2/pi).*(exp(-(a-1).^2.*snr/2)+exp(-(a+1).^2.*snr/2)- ... % -x
	  2*exp(-snr/2))-snr.*(1-2*Q(sqrt(snr))-Q((a-1).*sqrt(snr))+ ... 
	  Q((a+1).*sqrt(snr))) ...
	  ...
	  + 0.5*(snr.^2+snr).*(1-Q((a+1).*sqrt(snr))-Q((a-1).*sqrt(snr))) ... % + 1/2 x
	  + snr/2.*sqrt(snr/2/pi).*(-(a+1).*exp(-(a-1).^2.*snr/2) + ... 
     (1-a).*exp(-(a+1).^2.*snr/2)) ...
	  ...    
	  -(snr.^4/12+snr.^3/2+snr.^2/4).*(1-Q((a+1).*sqrt(snr)) - ... % - 1/12 x^4
     Q((a-1).*sqrt(snr))) + snr.^3/24.*sqrt(snr/2/pi).*(exp(-(a+1).^2 ...
	  /2.*snr).*(2*a.^3-2*a.^2+(6./snr+2).*a-(10./snr+2)) - ... 
     exp(-(a-1).^2/2.*snr).*(-2*a.^3-2*a.^2-(6./snr+2).*a-(10./snr+2))) ...
     ...
     +(snr.^6/45+snr.^5/3+snr.^4+snr.^3/3).*(1-Q((a+1).*sqrt(snr)) - ... % + 1/45 x^6
     Q((a-1).*sqrt(snr))) + sqrt(snr/2/pi)/45.*(exp(-(a+1).^2/2.*snr).* ...
        (snr.^5 + 14*snr.^4 + 33*snr.^3 - a.*(snr.^5 + 12*snr.^4 + 15*snr.^3) + ...
        a.^2.*(snr.^5 + 9*snr.^4) - a.^3.*(snr.^5 + 5*snr.^4) + (a.^4 - a.^5).*snr.^5) - ...
     exp(-(a-1).^2/2.*snr).*(snr.^5 + 14*snr.^4 + 33*snr.^3 + ...
        a.*(snr.^5 + 12*snr.^4 + 15*snr.^3) + a.^2.*(snr.^5 + 9*snr.^4) ...
        + a.^3.*(snr.^5 + 5*snr.^4) + (a.^4 + a.^5).*snr.^5)) ...
     ...
     -17/2520*(snr.^8+28*snr.^7+210*snr.^6+420*snr.^5+105*snr.^4).* ... % - 17/2520 x^8
     (1-Q((a+1).*sqrt(snr))-Q((a-1).*sqrt(snr))) - sqrt(snr/2/pi)*17/2520.* ...
      (exp(-(a+1).^2/2.*snr).*(snr.^7 + 27*snr.^6 + 185*snr.^5 + 279*snr.^4 ...
        - a.*(snr.^7 + 25*snr.^6 + 141*snr.^5 + 105*snr.^4) + a.^2.*(snr.^7 + 22*snr.^6 + 87*snr.^5) ...
        - a.^3.*(snr.^7 + 18*snr.^6 + 35*snr.^5) + a.^4.*(snr.^7 + 13*snr.^6) ...
        - a.^5.*(snr.^7 + 7*snr.^6) + (a.^6 - a.^7).*snr.^7) ...
     - exp(-(a-1).^2/2.*snr).*(snr.^7 + 27*snr.^6 + 185*snr.^5 + 279*snr.^4 ...
        + a.*(snr.^7 + 25*snr.^6 + 141*snr.^5 + 105*snr.^4) + a.^2.*(snr.^7 + 22*snr.^6 + 87*snr.^5) ...
        + a.^3.*(snr.^7 + 18*snr.^6 + 35*snr.^5) + a.^4.*(snr.^7 + 13*snr.^6) ...
        + a.^5.*(snr.^7 + 7*snr.^6) + (a.^6 + a.^7).*snr.^7)) ...
        ...
       +31/14175*(snr.^10+45*snr.^9+630*snr.^8+3150*snr.^7+4725*snr.^6+945*snr.^5).* ...% + 31/14175 x^10
         (1-Q((a+1).*sqrt(snr))-Q((a-1).*sqrt(snr))) + 31/14175*sqrt(snr/2/pi).* ...
       (exp(-(a+1).^2/2.*snr).*(snr.^9+44*snr.^8+588*snr.^7+2640*snr.^6+2895*snr.^5 ...
          - a.*(snr.^9+42*snr.^8+510*snr.^7+1830*snr.^6+945*snr.^5) ...
          + a.^2.*(snr.^9+39*snr.^8+405*snr.^7+975*snr.^6) - a.^3.*(snr.^9+35*snr.^8+285*snr.^7+315*snr.^6) ...
          + a.^4.*(snr.^9+30*snr.^8+165*snr.^7) - a.^5.*(snr.^9+24*snr.^8+63*snr.^7) ...
          + a.^6.*(snr.^9+17*snr.^8) - a.^7.*(snr.^9+9*snr.^8) ...
          + a.^8.*snr.^9 - a.^9.*snr.^9) ...
       - exp(-(a-1).^2/2.*snr).*(snr.^9+44*snr.^8+588*snr.^7+2640*snr.^6+2895*snr.^5 ...
          + a.*(snr.^9+42*snr.^8+510*snr.^7+1830*snr.^6+945*snr.^5) ...
          + a.^2.*(snr.^9+39*snr.^8+405*snr.^7+975*snr.^6) + a.^3.*(snr.^9+35*snr.^8+285*snr.^7+315*snr.^6) ...
          + a.^4.*(snr.^9+30*snr.^8+165*snr.^7) + a.^5.*(snr.^9+24*snr.^8+63*snr.^7) ...
          + a.^6.*(snr.^9+17*snr.^8) + a.^7.*(snr.^9+9*snr.^8) ...
          + a.^8.*snr.^9 + a.^9.*snr.^9));
d_taylow4 = -d_taylow4/log(2);

snr = hisnr;
a = 0.1832./snr;
d_taylow = - log(2) - snr.*(1-Q(-sqrt(snr)) + Q(sqrt(snr))) ...
	  + sqrt(2*snr/pi).*exp(-snr/2) + Q((a+1).*sqrt(snr)) + exp(4*snr) ...
	  .* Q((a+3).*sqrt(snr)) + ... % orignal bound
     ...
     log(2)*(1-Q((a-1).*sqrt(snr))-Q((a+1).*sqrt(snr))); % ln(2)
d_taylow = -d_taylow/log(2);

lower = [d_taylow4 d_taylow]-capacity;