WSJT-X/lib/fsk4hf/spb.m
Steven Franke db1454e927 Add reference for spb calculation.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8405 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2018-01-12 22:48:33 +00:00

90 lines
2.6 KiB
Mathematica

clear all;
global N
global R
global A
#-------------------------------------------------------------------------------
function retval = f1(theta)
global N;
global R;
retval=0.0;
gterm = gammaln(N/2) - gammaln((N+1)/2) - log(2*sqrt(pi));
rhs = -N*R*log(2);
lhs=gterm + (N-1)*log(sin(theta)) + log(1-(tan(theta).^2)/N) - log(cos(theta));
retval = rhs-real(lhs);
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function retval = d(N,i,x)
t1=(x.^2)/2;
t2=gammaln(N/2);
t3=-gammaln(i/2+1);
t4=-gammaln(N-i);
t5=(N-1-i)*log(sqrt(2)*x);
t6=-log(2)/2;
t7arg=1+(-1)^i * gammainc((x.^2)/2.0,(i+1)/2);
t7=log(t7arg);
retval=t1+t2+t3+t4+t5+t6+t7;
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function retval = maxstar(x1,x2)
retval = max(x1,x2)+log(1+exp(-abs(x1-x2)));
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function retval = spb_integrand(x)
global N;
global A;
t1=log(N-1);
t2=-N*(A^2)/2;
t3=-0.5*log(2*pi);
t4=(N-2)*log(sin(x));
arg=sqrt(N)*A*cos(x);
t5=maxstar(d(N,0,arg),d(N,1,arg));
for i=2:N-1
t5=maxstar(t5,d(N,i,arg));
endfor
retval=exp(t1+t2+t3+t4+t5);
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function retval = qfunc(x)
retval = 0.5 * erfc(x/sqrt(2));
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# Calculate sphere packing lower bound on the probability of word error
# given block length (N), code rate (R), and Eb/No.
#
# Ref:
# "Log-Domain Calculation of the 1959 Sphere-Packing Bound with Application to
# M-ary PSK Block Coded Modulation," Igal Sason and Gil Weichman,
# doi: 10.1109/EEEI.2006.321097
#-------------------------------------------------------------------------------
N=174
K=75
R=K/N
delta=0.01;
[ths,fval,info,output]=fzero(@f1,[delta,pi/2-delta], optimset ("jacobian", "off"));
for ebnodb=-6:0.5:4
ebno=10^(ebnodb/10.0);
esno=ebno*R;
A=sqrt(2*esno);
term1=quadcc(@spb_integrand,ths,pi/2);
term2=qfunc(sqrt(N)*A);
pe=term1+term2;
ps=1-pe;
printf("%f %f\n",ebnodb,ps);
endfor