w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2024-11-22 04:11:16 -05:00
Code Issues Packages Projects Releases Wiki Activity

Add spb.m - a gnu octave script to calculate the sphere packing bound. Only tested for low SNR and short block lengths.

Browse Source 
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8399 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Steven Franke 2018-01-10 19:43:34 +00:00
parent 4f0e73f68f
commit 850518aa27
1 changed files with 84 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

84
lib/fsk4hf/spb.m Normal file
View File

@ -0,0 +1,84 @@
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.
#-------------------------------------------------------------------------------
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