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;