simulation to plot error curves for BIll VK5DSP's short LDPC codes

git-svn-id: https://svn.code.sf.net/p/freetel/code@3059 01035d8c-6547-0410-b346-abe4f91aad63

diff --git a/codec2-dev/octave/ldpc_short.m b/codec2-dev/octave/ldpc_short.m
new file mode 100644 (file)
index 0000000..22559e1
--- /dev/null
+++ b/codec2-dev/octave/ldpc_short.m
@@ -0,0 +1,142 @@
+% ldpc.m\r
+%\r
+% David Rowe Mar 2017\r
+% Based on Bill Cowley's LDPC simulations\r
+%\r
+% Octave simulation of BPSK with short LDPC codes developed by Bill.  First step\r
+% in use of LDPC codes with FreeDV and Codec 2 700C.\r
+%\r
+% See lpdc.m for instruction son how to install CML library\r
+\r
+1;\r
+\r
+function init_cml\r
+  currentdir = pwd;\r
+  thiscomp = computer;\r
+\r
+  if strfind(thiscomp, 'pc-linux-gnu')==8 \r
+     if exist('LdpcEncode')==0, \r
+        %cd '~/cml'\r
+        cd '/home/david/Desktop/cml'\r
+        CmlStartup      \r
+     end\r
+  end\r
+  cd(currentdir); \r
+end\r
+\r
+function sim_out = ldpc4(HRA, sim_in, genie_Es);\r
+\r
+  estEsN0 = 0\r
+\r
+  framesize = sim_in.framesize;\r
+  rate      = sim_in.rate;\r
+  mod_order = sim_in.mod_order;\r
+  Ntrials   = sim_in.Ntrials;\r
+  Esvec     = sim_in.Esvec;\r
+\r
+  demod_type = 0;\r
+  decoder_type = 0;\r
+  max_iterations = 100;\r
+  bps = code_param.bits_per_symbol = log2(mod_order);\r
+\r
+  [H_rows, H_cols] = Mat2Hrows(HRA); \r
+  code_param.H_rows = H_rows; \r
+  code_param.H_cols = H_cols;\r
+  code_param.P_matrix = [];\r
+  code_param.data_bits_per_frame = length(code_param.H_cols) - length( code_param.P_matrix ); \r
+\r
+\r
+  for ne = 1:length(Esvec)\r
+    Es = Esvec(ne);\r
+    EsNo = 10^(Es/10);\r
+    EbNodB = Es - 10*log10(code_param.bits_per_symbol * rate);\r
+\r
+    Terrs = 0;  Tbits = 0;  Ferrs = 0;\r
+\r
+    for nn = 1: Ntrials        \r
+      data = round( rand( 1, code_param.data_bits_per_frame ) );\r
+      codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix );\r
+      code_param.code_bits_per_frame = length( codeword );\r
+      Nsymb = code_param.code_bits_per_frame/bps;      \r
+       \r
+      % modulate\r
+      s = 1 - 2 * codeword;   \r
+      code_param.symbols_per_frame = length( s );\r
+              \r
+      variance = 1/(2*EsNo);\r
+      noise = sqrt(variance)* randn(1,code_param.symbols_per_frame); \r
+      r = s + noise;\r
+\r
+      Nr = length(r);  \r
+      % in the binary case the LLRs are just a scaled version of the rx samples ...\r
+      % if the EsNo is known -- use the following \r
+      if (genie_Es) \r
+       input_decoder_c = 4 * EsNo * r;   \r
+      else\r
+        r = r / mean(abs(r));       % scale for signal unity signal  \r
+       estvar = var(r-sign(r)); \r
+       estEsN0 = 1/(2* estvar); \r
+       input_decoder_c = 4 * estEsN0 * r;\r
+      end\r
+\r
+      [x_hat, PCcnt] = MpDecode( input_decoder_c, code_param.H_rows, code_param.H_cols, ...\r
+                                 max_iterations, decoder_type, 1, 1);\r
+      Niters = sum(PCcnt!=0);\r
+      detected_data = x_hat(Niters,:);\r
+      error_positions = xor( detected_data(1:code_param.data_bits_per_frame), data );\r
+      Nerrs = sum( error_positions);\r
+        \r
+      if Nerrs>0,  Ferrs = Ferrs +1;  end\r
+      Terrs = Terrs + Nerrs;\r
+      Tbits = Tbits + code_param.data_bits_per_frame;\r
+    end\r
+      \r
+    printf("EbNo: %3.1f dB BER: %5.4f Nbits: %4d Nerrs: %4d\n", EbNodB, Terrs/Tbits, Tbits, Terrs);\r
+\r
+    TERvec(ne) = Terrs;\r
+    FERvec(ne) = Ferrs;\r
+    BERvec(ne) = Terrs/ Tbits;\r
+    Ebvec = Esvec - 10*log10(code_param.bits_per_symbol * rate);\r
+    \r
+    sim_out.BERvec = BERvec;\r
+    sim_out.Ebvec = Ebvec;\r
+    sim_out.FERvec = FERvec;\r
+    sim_out.TERvec  = TERvec;\r
+    sim_out.cpumins = cputime/60;    \r
+  end\r
+endfunction\r
+\r
+\r
+more off;\r
+format;\r
+init_cml\r
+\r
+sim_in.Esvec = -3:0.5:6; \r
+load Hs_112_112.mat\r
+\r
+HRA = H2;\r
+\r
+[Nr Nc] = size(HRA);  \r
+\r
+sim_in.rate = (Nc-Nr)/Nc\r
+sim_in.framesize = Nc;\r
+\r
+sim_in.mod_order = 2; \r
+sim_in.modulation = 'BPSK';\r
+sim_in.mapping = 'gray';\r
+\r
+sim_in.Lim_Ferrs= 100;\r
+sim_in.Ntrials =  500;\r
+\r
+sim_out = ldpc4(HRA, sim_in, 1);\r
+\r
+figure(1)\r
+semilogy(sim_out.Ebvec,  sim_out.BERvec)\r
+xlabel('Eb/N0')\r
+ylabel('BER')\r
+grid\r
+figure(2)\r
+semilogy(sim_out.Ebvec,  sim_out.FERvec/sim_in.Ntrials ,  col)\r
+xlabel('Eb/N0')\r
+ylabel('PER')\r
+grid\r