llr calculations in C, modified Octave to use bit exact procedure

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

diff --git a/codec2-dev/octave/ldpc_fsk_lib.m b/codec2-dev/octave/ldpc_fsk_lib.m
index c9df72103c43d0c4aa56c8c0e4e91d0c8a603e2b..663e9d94720976ce215ac03fdbb27cfebad6368d 100644 (file)
--- a/codec2-dev/octave/ldpc_fsk_lib.m
+++ b/codec2-dev/octave/ldpc_fsk_lib.m
@@ -141,23 +141,44 @@ function codeword = ldpc_encode(code_param, data)
 endfunction\r
 \r
 \r
+function llr = sd_to_llr(sd)\r
+    %printf("mean: %2.17g\n", mean(abs(sd)));\r
+    sd = sd / mean(abs(sd));\r
+    x = sd - sign(sd);\r
+    sumsq = sum(x.^2);\r
+    summ = sum(x);\r
+    mn = summ/length(sd);\r
+    estvar = sumsq/length(sd) - mn*mn; \r
+    %printf("estvar: %2.17g\n", estvar);\r
+    estEsN0 = 1/(2* estvar + 1E-3); \r
+    llr = 4 * estEsN0 * sd;\r
+    %printf("\n");\r
+    %for i=st:st+4\r
+    %  printf("%2.17g ", input_decoder_c(i))\r
+    %end\r
+    %printf("\n");\r
+endfunction\r
+\r
 % LDPC decoder - note it estimates EsNo from received symbols\r
 \r
 function [detected_data Niters] = ldpc_decode(r, code_param, max_iterations, decoder_type)\r
   % in the binary case the LLRs are just a scaled version of the rx samples ..\r
 \r
+ #{\r
   r = r / mean(abs(r));       % scale for signal unity signal  \r
   estvar = var(r-sign(r)); \r
   estEsN0 = 1/(2* estvar + 1E-3); \r
   input_decoder_c = 4 * estEsN0 * r;\r
+ #}\r
+  llr = sd_to_llr(r);\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
+  [x_hat, PCcnt] = MpDecode(llr, 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
-\r
+  \r
   if isfield(code_param, "c_include_file")\r
-    write_code_to_C_include_file(code_param, max_iterations, decoder_type, input_decoder_c, x_hat);\r
+    write_code_to_C_include_file(code_param, max_iterations, decoder_type, llr, x_hat);\r
   end\r
 end\r
 \r
@@ -210,7 +231,7 @@ function write_code_to_C_include_file(code_param, max_iterations, decoder_type,
 \r
   fprintf(f,"\ndouble input[] = {\n");\r
   for i=1:length(input_decoder_c)\r
-    fprintf(f, "%f", input_decoder_c(i));\r
+    fprintf(f, "%.17g", input_decoder_c(i));\r
     if i == length(input_decoder_c)\r
       fprintf(f,"\n};\n");\r
     else\r

diff --git a/codec2-dev/octave/test_ldpc_fsk_lib.m b/codec2-dev/octave/test_ldpc_fsk_lib.m
index f86a8b936eb56d4eb5e3ffdc51dbafd4b195342f..e7be8b9d7469f9418778d93cae1069712e6a4c30 100644 (file)
--- a/codec2-dev/octave/test_ldpc_fsk_lib.m
+++ b/codec2-dev/octave/test_ldpc_fsk_lib.m
@@ -14,7 +14,7 @@ function [data code_param] = simple_ut(c_include_file)
   HRA = full(HRA);  \r
   max_iterations = 100;\r
   decoder_type = 0;\r
-  EsNodB = 4;\r
+  EsNodB = 3;\r
   mod_order = 2;\r
 \r
   code_param = ldpc_init(HRA, mod_order);\r
@@ -152,6 +152,87 @@ function test_c_encoder
 end\r
 \r
 \r
+function test_c_decoder\r
+  load('H2064_516_sparse.mat');\r
+  HRA = full(HRA);  \r
+  max_iterations = 100;\r
+  decoder_type = 0;\r
+  mod_order = 2;\r
+  frames = 10;\r
+  EsNodB = 2;\r
+  sdinput = 1;\r
+\r
+  EsNo = 10^(EsNodB/10);\r
+  variance = 1/(2*EsNo);\r
+\r
+  code_param = ldpc_init(HRA, mod_order);\r
+  data = round(rand(1,code_param.data_bits_per_frame*frames));\r
+  \r
+  f = fopen("data.bin","wt"); fwrite(f, data, "uint8"); fclose(f);\r
+  system("../src/ldpc_enc data.bin codewords.bin"); \r
+  f = fopen("codewords.bin","rb"); codewords = fread(f, "uint8")'; fclose(f);\r
+\r
+  s = 1 - 2 * codewords;   \r
+  noise = sqrt(variance)*randn(1,code_param.symbols_per_frame*frames); \r
+  r = s + noise;\r
+\r
+  % calc LLRs frame by frame\r
+\r
+  for i=1:frames\r
+    st = (i-1)*code_param.symbols_per_frame+1;\r
+    en = st + code_param.symbols_per_frame-1;\r
+    llr(st:en) = sd_to_llr(r(st:en));    \r
+  end\r
+\r
+  % Outboard C decoder\r
+\r
+  if sdinput\r
+    f = fopen("sd.bin","wb"); fwrite(f, r, "double"); fclose(f);\r
+    system("../src/ldpc_dec sd.bin data_out.bin --sdinput"); \r
+  else\r
+    f = fopen("llr.bin","wb"); fwrite(f, llr, "double"); fclose(f);\r
+    system("../src/ldpc_dec llr.bin data_out.bin"); \r
+  end\r
+\r
+  f = fopen("data_out.bin","rb"); data_out = fread(f, "uint8")'; fclose(f);\r
+  \r
+  Nerrs = Nerrs2 = zeros(1,frames);;\r
+  for i=1:frames\r
+\r
+    % Check C decoder\r
+    \r
+    data_st = (i-1)*code_param.data_bits_per_frame+1;\r
+    data_en = data_st+code_param.data_bits_per_frame-1;\r
+    st = (i-1)*code_param.symbols_per_frame+1;\r
+    en = st+code_param.data_bits_per_frame-1;\r
+    data_out_c = data_out(st:en);\r
+    error_positions = xor(data_out_c, data(data_st:data_en));\r
+    Nerrs(i) = sum(error_positions);\r
+\r
+    % Octave decoder \r
+\r
+    st = (i-1)*code_param.symbols_per_frame+1; en = st+code_param.symbols_per_frame-1;\r
+    [detected_data Niters] = ldpc_decode(r(st:en), code_param, max_iterations, decoder_type);\r
+    st = (i-1)*code_param.data_bits_per_frame+1; en = st+code_param.data_bits_per_frame-1;\r
+    data_out_octave = detected_data(1:code_param.data_bits_per_frame);\r
+    error_positions = xor(data_out_octave, data(st:en));\r
+    Nerrs2(i) = sum(error_positions);\r
+    %printf("%4d ", Niters);\r
+  end\r
+  printf("Errors per frame:\nC.....:");\r
+  for i=1:frames\r
+    printf("%4d ", Nerrs(i));\r
+  end\r
+  printf("\nOctave:");\r
+  for i=1:frames\r
+    printf("%4d ", Nerrs2(i));\r
+  end\r
+  printf("\n");\r
+\r
+end\r
+\r
+\r
+\r
 % Start simulation --------------------------------------------------------\r
 \r
 more off;\r
@@ -185,7 +266,7 @@ rand('state',1);
 \r
 % binary flags to run various demos, e.g. "15" to run 1 .. 8\r
 \r
-demos = 32;\r
+demos = 64;\r
 \r
 if bitand(demos,1)\r
   printf("simple_ut....\n");\r
@@ -230,7 +311,10 @@ if bitand(demos,16)
   printf("Nerrs = %d\n", Nerrs);\r
 end\r
 \r
-\r
 if bitand(demos,32)\r
   test_c_encoder;\r
 end\r
+\r
+if bitand(demos,64)\r
+  test_c_decoder;\r
+end\r