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;