first pass a AWGN perf curves for Bill's short codes

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

diff --git a/codec2-dev/octave/HRA_112_112.txt b/codec2-dev/octave/HRA_112_112.txt
index e446c4a9e1bf9e968eb82fb0c788aa9cc0d8b0f5..e7746401261a165f1959a52e0c8afa9f8c4d1240 100644 (file)
--- a/codec2-dev/octave/HRA_112_112.txt
+++ b/codec2-dev/octave/HRA_112_112.txt
@@ -1,5 +1,5 @@
 # Created by Octave 3.8.1, Wed Mar 15 17:24:05 2017 ACDT <david@penetrator>
-# name: H2
+# name: HRA_112_112
 # type: matrix
 # rows: 112
 # columns: 224

diff --git a/codec2-dev/octave/HRA_112_56.txt b/codec2-dev/octave/HRA_112_56.txt
index e22ce658dbf73f5cf000699aa9408917452b6f6a..1df20d8cdc2774d854876639fc15e5be4a441375 100644 (file)
--- a/codec2-dev/octave/HRA_112_56.txt
+++ b/codec2-dev/octave/HRA_112_56.txt
@@ -1,5 +1,5 @@
 # Created by Octave 3.8.1, Wed Mar 15 17:59:13 2017 ACDT <david@penetrator>
-# name: H1
+# name: HRA_112_56
 # type: matrix
 # rows: 56
 # columns: 168

diff --git a/codec2-dev/octave/HRA_56_28.txt b/codec2-dev/octave/HRA_56_28.txt
index f31c00faebe83467b53b3a9e9ca77a9cf98a72b2..ae7624849f99a979992ce58e72efb9850df8db7b 100644 (file)
--- a/codec2-dev/octave/HRA_56_28.txt
+++ b/codec2-dev/octave/HRA_56_28.txt
@@ -1,5 +1,5 @@
 # Created by Octave 3.8.1, Wed Mar 15 18:20:20 2017 ACDT <david@penetrator>
-# name: H1
+# name: HRA_56_28
 # type: matrix
 # rows: 28
 # columns: 84

diff --git a/codec2-dev/octave/HRA_56_56.txt b/codec2-dev/octave/HRA_56_56.txt
index b35e0e4d48dd058337725e7a7527a412efdbac50..8a18f2d611fe874dc3b94c4f5a898e952ab41752 100644 (file)
--- a/codec2-dev/octave/HRA_56_56.txt
+++ b/codec2-dev/octave/HRA_56_56.txt
@@ -1,5 +1,5 @@
 # Created by Octave 3.8.1, Wed Mar 15 18:12:09 2017 ACDT <david@penetrator>
-# name: H1
+# name: HRA_56_56
 # type: matrix
 # rows: 56
 # columns: 112

diff --git a/codec2-dev/octave/ldpc_short.m b/codec2-dev/octave/ldpc_short.m
index e769bf740dec72d540f1c6bcb892e605c0944e2a..a60924d35999ca6cdfb45a1dd38a0c1f49338699 100644 (file)
--- a/codec2-dev/octave/ldpc_short.m
+++ b/codec2-dev/octave/ldpc_short.m
@@ -25,36 +25,50 @@ function init_cml
 end\r
 \r
 \r
-function sim_out = ldpc4(HRA, Ntrials, Esvec, genie_Es);\r
-\r
-  [Nr Nc] = size(HRA);  \r
+function sim_out = run_sim(HRA, Ntrials, Esvec, genie_Es, packet_size, golay = 0)\r
+\r
+  if golay\r
+    rate = 0.5;\r
+    code_param.data_bits_per_frame = 12;\r
+    framesize = 24;\r
+  else\r
+    % LDPC code\r
+\r
+    [Nr Nc] = size(HRA);  \r
+    rate = (Nc-Nr)/Nc;\r
+    framesize = Nc;\r
+    modulation = 'BPSK';\r
+    demod_type = 0;\r
+    decoder_type = 0;\r
+    max_iterations = 100;\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
+  end\r
 \r
-  rate = (Nc-Nr)/Nc;\r
-  framesize = Nc;\r
   mod_order = 2; \r
-  modulation = 'BPSK';\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
+  % loop around each Esvec point\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
+    Terrs = 0;  Tbits = 0;  Ferrs = 0; \r
+    tx_bits = rx_bits = [];\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
+      tx_bits = [tx_bits data];\r
+      if golay\r
+        codeword = egolayenc(data);\r
+      else\r
+        codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix );\r
+      end\r
       code_param.code_bits_per_frame = length( codeword );\r
       Nsymb = code_param.code_bits_per_frame/bps;      \r
        \r
@@ -67,37 +81,64 @@ function sim_out = ldpc4(HRA, Ntrials, Esvec, genie_Es);
       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
+\r
+      if golay\r
+        detected_data = egolaydec(r < 0);\r
+        detected_data = detected_data(code_param.data_bits_per_frame+1:framesize);\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
+\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
+\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
+        detected_data = detected_data(1:code_param.data_bits_per_frame);\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
+      rx_bits = [rx_bits detected_data];\r
+\r
+      error_positions = xor( detected_data, 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
+\r
     end\r
       \r
-    printf("EbNo: %3.1f dB BER: %5.4f Nbits: %4d Nerrs: %4d\n", EbNodB, Terrs/Tbits, Tbits, Terrs);\r
+    % count packet errors using supplied packet size.  Operating one\r
+    % one big string of bits lets us account for cases were FEC framesize\r
+    % is less than packet size\r
+\r
+    Perrs = 0; Tpackets = 0;\r
+    for i=1:packet_size:length(tx_bits)-packet_size\r
+      error_positions = xor( tx_bits(i:i+packet_size-1), rx_bits(i:i+packet_size-1));\r
+      Nerrs = sum(error_positions);\r
+      if Nerrs>0,  Perrs = Perrs +1;  end\r
+      Tpackets++;\r
+    end\r
+\r
+    printf("EbNo: %3.1f dB BER: %5.4f PER: %5.4f Nbits: %4d Nerrs: %4d Tpackets: %4d Perr: %4d\n", EbNodB, Terrs/Tbits, Perrs/ Tpackets, Tbits, Terrs, Tpackets, Perrs);\r
 \r
     TERvec(ne) = Terrs;\r
     FERvec(ne) = Ferrs;\r
     BERvec(ne) = Terrs/ Tbits;\r
+    PERvec(ne) = Perrs/ Tpackets;\r
     Ebvec = Esvec - 10*log10(code_param.bits_per_symbol * rate);\r
     \r
     sim_out.BERvec = BERvec;\r
+    sim_out.PERvec = PERvec;\r
     sim_out.Ebvec = Ebvec;\r
     sim_out.FERvec = FERvec;\r
     sim_out.TERvec  = TERvec;\r
@@ -107,42 +148,67 @@ endfunction
 \r
 % Start simulation here ----------------------------------------------\r
 \r
+rand('seed',1);\r
+randn('seed',1);\r
 more off;\r
 format;\r
-init_cml\r
+init_cml;\r
 \r
-Ntrials =  2000;\r
+Ntrials = 500;\r
 Esvec = -3:0.5:3; \r
+packet_size = 28;\r
+\r
+\r
+load HRA_112_112.txt\r
+load HRA_112_56.txt\r
+load HRA_56_56.txt\r
+load HRA_56_28.txt\r
+\r
+sim_out1 = run_sim(HRA_112_112, Ntrials, Esvec, 1, packet_size);\r
+sim_out2 = run_sim(HRA_112_56 , Ntrials, Esvec, 1, packet_size);\r
+sim_out3 = run_sim(HRA_56_56  , Ntrials*2, Esvec, 1, packet_size);\r
+sim_out4 = run_sim(HRA_56_28  , Ntrials*2, Esvec, 1, packet_size);\r
+sim_out5 = run_sim([], Ntrials*10, Esvec, 1, packet_size, 1);\r
 \r
-sim_out1 = ldpc4(H1, Ntrials, Esvec, 1);\r
-sim_out2 = ldpc4(H2, Ntrials, Esvec, 1);\r
-sim_out3 = ldpc4(H3, Ntrials, Esvec, 1);\r
-sim_out4 = ldpc4(H4, Ntrials, Esvec, 1);\r
+Ebvec_theory = -2:0.5:6;\r
+uncoded_BER_theory = 0.5*erfc(sqrt(10.^(Ebvec_theory/10)));\r
 \r
-EbNodBvec = sim_out1.Ebvec;\r
-uncoded_BER_theory = 0.5*erfc(sqrt(10.^(EbNodBvec/10)));\r
-uncoded_PER_theory = uncoded_BER_theory*sim_out1.framesize;\r
+% need standard packet size to compare\r
+% packet error if bit 0, or bit 1, or bit 2 .....\r
+%              or bit 0 and bit 1\r
+% no packet error if all bits ok (1-p(0))*(1-p(1))\r
+% P(packet error) = p(0)+p(1)+....\r
+\r
+uncoded_PER_theory = 1 - (1-uncoded_BER_theory).^packet_size;\r
 \r
 figure(1); clf;\r
-semilogy(EbNodBvec,  uncoded_BER_theory, 'b')\r
+semilogy(Ebvec_theory,  uncoded_BER_theory, 'b+-;BPSK theory;','markersize', 10, 'linewidth', 2)\r
 hold on;\r
-semilogy(EbNodBvec,  sim_out1.BERvec, 'g')\r
-semilogy(EbNodBvec,  sim_out2.BERvec, 'r')\r
-semilogy(EbNodBvec,  sim_out3.BERvec, 'c')\r
-semilogy(EbNodBvec,  sim_out4.BERvec, 'k')\r
+semilogy(sim_out1.Ebvec, sim_out1.BERvec, 'g+-;rate 1/2 HRA 112 112;','markersize', 10, 'linewidth', 2)\r
+semilogy(sim_out2.Ebvec, sim_out2.BERvec, 'r+-;rate 2/3 HRA 112 56;','markersize', 10, 'linewidth', 2)\r
+semilogy(sim_out3.Ebvec, sim_out3.BERvec, 'c+-;rate 1/2 HRA 56 56;','markersize', 10, 'linewidth', 2)\r
+semilogy(sim_out4.Ebvec, sim_out4.BERvec, 'k+-;rate 2/3 HRA 56 28;','markersize', 10, 'linewidth', 2)\r
+semilogy(sim_out5.Ebvec, sim_out5.BERvec, 'm+-;rate 1/2 Golay (24,12);','markersize', 10, 'linewidth', 2)\r
 hold off;\r
 xlabel('Eb/No')\r
 ylabel('BER')\r
 grid\r
+legend("boxoff");\r
+\r
 \r
 figure(2); clf;\r
-semilogy(EbNodBvec,  uncoded_PER_theory, 'b')\r
+semilogy(Ebvec_theory,  uncoded_PER_theory, 'b+-;BPSK theory;','markersize', 10, 'linewidth', 2)\r
 hold on;\r
-semilogy(EbNodBvec,  sim_out1.FERvec/Ntrials, 'g')\r
-semilogy(EbNodBvec,  sim_out2.FERvec/Ntrials, 'r')\r
-semilogy(EbNodBvec,  sim_out3.FERvec/Ntrials, 'c')\r
-semilogy(EbNodBvec,  sim_out4.FERvec/Ntrials, 'k')\r
+semilogy(sim_out1.Ebvec, sim_out1.PERvec, 'g+-;rate 1/2 HRA 112 112;','markersize', 10, 'linewidth', 2)\r
+semilogy(sim_out2.Ebvec, sim_out2.PERvec, 'r+-;rate 2/3 HRA 112 56;','markersize', 10, 'linewidth', 2)\r
+semilogy(sim_out3.Ebvec, sim_out3.PERvec, 'c+-;rate 1/2 HRA 56 56;','markersize', 10, 'linewidth', 2)\r
+semilogy(sim_out4.Ebvec, sim_out4.PERvec, 'k+-;rate 2/3 HRA 56 28;','markersize', 10, 'linewidth', 2)\r
+semilogy(sim_out5.Ebvec, sim_out5.PERvec, 'm+-;rate 1/2 Golay (24,12);','markersize', 10, 'linewidth', 2)\r
 hold off;\r
 xlabel('Eb/No')\r
 ylabel('PER')\r
 grid\r
+legend("boxoff");\r
+\r
+\r
+\r