function [upper] = BINCAPUP(SNR, capacity)
% BINCAPUP - Upper bound on binary input AWGN channel capacity
%   C = BINCAPUP(SNR) calculates an upper 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.305./snr;
cap_high = - log(2) - 2*snr.*Q(sqrt(snr)) + sqrt(2*snr/pi).*exp(-snr/2)...
     + log(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^2
	  + 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));
cap_high = -cap_high/log(2);
  
snr = hisnr;
a = 0.3465./snr;
d_tayhi = - log(2) - snr.*(1-Q(-sqrt(snr)) + Q(sqrt(snr))) ...
	  + sqrt(2*snr/pi).*exp(-snr/2) + ...
	  log(2)*(Q((a+1).*sqrt(snr)) + exp(4*snr).*Q((a+3).*sqrt(snr))) + ... %orig.
	  ...
     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)));
d_tayhi = -d_tayhi/log(2);

upper = [cap_high d_tayhi]-capacity;