From 91c3386e1c3e5c1b2301997c68a0b4e23749bcdb Mon Sep 17 00:00:00 2001
From: baobrien <baobrien@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Sun, 27 Sep 2015 02:42:58 +0000
Subject: [PATCH] fsk4 ber curves work now

git-svn-id: https://svn.code.sf.net/p/freetel/code@2404 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/octave/fsk4.m | 43 +++++++++++++++++++++++++++++++++-------
 1 file changed, 36 insertions(+), 7 deletions(-)

diff --git a/codec2-dev/octave/fsk4.m b/codec2-dev/octave/fsk4.m
index a391406f..8da6cfd2 100644
--- a/codec2-dev/octave/fsk4.m
+++ b/codec2-dev/octave/fsk4.m
@@ -4,7 +4,7 @@
 graphics_toolkit("gnuplot");
 
 fm;
-
+pkg load signal;
 % Frequency response of the DMR raised cosine filter 
 % from ETSI TS 102 361-1 V2.2.1 page 111
 dmr.tx_filt_resp = @(f) 1.0*(f<=1920) - cos((pi*f)/1920).*1.0.*(f>1920 & f<=2880);
@@ -232,13 +232,13 @@ endfunction
 % Bit error rate test
 % for a noise-free channel
 % now supports noisy channels
-function [ber thrber] = nfbert(aEsNodB)
+function [ber thrcoh thrncoh] = nfbert(aEsNodB)
   global dmr_info;
   global nxdn_info;
   rand('state',1); 
   randn('state',1);
 
-  bitcnt = 96000;
+  bitcnt = 480000;
   test_bits = [zeros(1,100) rand(1,bitcnt)>.5]; %Random bits. Pad with zeros to prime the filters
   fsk4_states.M = 1;
   fsk4_states = fsk4_init(fsk4_states,2400,dmr_info);
@@ -273,23 +273,52 @@ function [ber thrber] = nfbert(aEsNodB)
   offset = ox;
   printf("\ncoarse timing: %d nerr: %d\n", offset, best_nerr);
 
-  thrber = erfc(sqrt(EbNo));
+  % Coherent BER theory
+  thrcoh = erfc(sqrt(EbNo));
+
+  % non-coherent BER theory calculation
+  % It was complicated, so I broke it up
+
+  ms = 4;
+  ns = (1:ms-1);
+  as = (-1).^(ns+1);
+  bs = (as./(ns+1));
+  
+  cs = ((ms-1)./ns);
+
+  ds = ns.*log2(ms);
+  es = ns+1;
+  fs = exp( -(ds./es)*EbNo );
+  
+  thrncoh = ((ms/2)/(ms-1)) * sum(bs.*((ms-1)./ns).*exp( -(ds./es)*EbNo ));
+  %thrncoh = (2/3)*sum( (((-1).^(ns+1))./(ns+1)) .* (3./ns) .* exp( -((ns.*2)./ns)*EbNo )
   %plot(xor(rx_bits(offset:length(rx_bits)),test_bits(1:length(rx_bits)+1-offset)));
   %plot((1:1000),rx_bits(1:1000),(1:1000),rx_err(1:1000));
 endfunction
 
-pkg load parallel
+
 function fsk4_ber_curves
   EbNodB = 4:13;
-  bers_tnco = bers_real = ones(1,length(EbNodB));
+  bers_tco = bers_real = bers_tnco = ones(1,length(EbNodB));
   %for ii=(1:length(EbNodB));
   %  [bers_real(ii),bers_tnco(ii)] = nfbert(EbNodB(ii));
   %end
-  [bers_real,bers_tnco] = pararrayfun(5,@nfbert,EbNodB);
+  try
+    pkg load parallel
+    [bers_real,bers_tco,bers_tnco] = pararrayfun(floor(1.25*nproc()),@nfbert,EbNodB);
+  catch
+    printf("You should install package parallel. It'll make this run way faster\n");
+    for ii=(1:length(EbNodB));
+      [bers_real(ii),bers,tco(ii),bers_tnco(ii)] = nfbert(EbNodB(ii));
+    end
+  end_try_catch
   figure;
   clf;
+  bers_tnco
   semilogy(EbNodB, bers_tnco,'r;4FSK non-coherent theory;')
   hold on;
+
+  semilogy(EbNodB, bers_tco,'b;4FSK coherent theory;')
   semilogy(EbNodB, bers_real ,'g;4FSK non-coherent sim;')
   hold off;
   grid("minor");
-- 
2.25.1