diff --git a/lib/fsk4hf/spb.m b/lib/fsk4hf/spb.m
new file mode 100644
index 000000000..cc8103aac
--- /dev/null
+++ b/lib/fsk4hf/spb.m
@@ -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