From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Wed, 14 Sep 2016 10:00:45 +0000 (+0000)
Subject: llr calculations in C, modified Octave to use bit exact procedure
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=07aea06e6d286ff7996a9452cd9dbca1f24030b6;p=freetel-svn-tracking.git

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 c9df7210..663e9d94 100644
--- 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
 
 
+function llr = sd_to_llr(sd)
+    %printf("mean: %2.17g\n", mean(abs(sd)));
+    sd = sd / mean(abs(sd));
+    x = sd - sign(sd);
+    sumsq = sum(x.^2);
+    summ = sum(x);
+    mn = summ/length(sd);
+    estvar = sumsq/length(sd) - mn*mn; 
+    %printf("estvar: %2.17g\n", estvar);
+    estEsN0 = 1/(2* estvar + 1E-3); 
+    llr = 4 * estEsN0 * sd;
+    %printf("\n");
+    %for i=st:st+4
+    %  printf("%2.17g ", input_decoder_c(i))
+    %end
+    %printf("\n");
+endfunction
+
 % LDPC decoder - note it estimates EsNo from received symbols
 
 function [detected_data Niters] = ldpc_decode(r, code_param, max_iterations, decoder_type)
   % in the binary case the LLRs are just a scaled version of the rx samples ..
 
+ #{
   r = r / mean(abs(r));       % scale for signal unity signal  
   estvar = var(r-sign(r)); 
   estEsN0 = 1/(2* estvar + 1E-3); 
   input_decoder_c = 4 * estEsN0 * r;
+ #}
+  llr = sd_to_llr(r);
 
-  [x_hat, PCcnt] = MpDecode( input_decoder_c, code_param.H_rows, code_param.H_cols, ...
-                             max_iterations, decoder_type, 1, 1);         
+  [x_hat, PCcnt] = MpDecode(llr, code_param.H_rows, code_param.H_cols, ...
+                            max_iterations, decoder_type, 1, 1);         
   Niters = sum(PCcnt!=0);
   detected_data = x_hat(Niters,:);
-
+  
   if isfield(code_param, "c_include_file")
-    write_code_to_C_include_file(code_param, max_iterations, decoder_type, input_decoder_c, x_hat);
+    write_code_to_C_include_file(code_param, max_iterations, decoder_type, llr, x_hat);
   end
 end
 
@@ -210,7 +231,7 @@ function write_code_to_C_include_file(code_param, max_iterations, decoder_type,
 
   fprintf(f,"\ndouble input[] = {\n");
   for i=1:length(input_decoder_c)
-    fprintf(f, "%f", input_decoder_c(i));
+    fprintf(f, "%.17g", input_decoder_c(i));
     if i == length(input_decoder_c)
       fprintf(f,"\n};\n");
     else
diff --git a/codec2-dev/octave/test_ldpc_fsk_lib.m b/codec2-dev/octave/test_ldpc_fsk_lib.m
index f86a8b93..e7be8b9d 100644
--- 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);  
   max_iterations = 100;
   decoder_type = 0;
-  EsNodB = 4;
+  EsNodB = 3;
   mod_order = 2;
 
   code_param = ldpc_init(HRA, mod_order);
@@ -152,6 +152,87 @@ function test_c_encoder
 end
 
 
+function test_c_decoder
+  load('H2064_516_sparse.mat');
+  HRA = full(HRA);  
+  max_iterations = 100;
+  decoder_type = 0;
+  mod_order = 2;
+  frames = 10;
+  EsNodB = 2;
+  sdinput = 1;
+
+  EsNo = 10^(EsNodB/10);
+  variance = 1/(2*EsNo);
+
+  code_param = ldpc_init(HRA, mod_order);
+  data = round(rand(1,code_param.data_bits_per_frame*frames));
+  
+  f = fopen("data.bin","wt"); fwrite(f, data, "uint8"); fclose(f);
+  system("../src/ldpc_enc data.bin codewords.bin"); 
+  f = fopen("codewords.bin","rb"); codewords = fread(f, "uint8")'; fclose(f);
+
+  s = 1 - 2 * codewords;   
+  noise = sqrt(variance)*randn(1,code_param.symbols_per_frame*frames); 
+  r = s + noise;
+
+  % calc LLRs frame by frame
+
+  for i=1:frames
+    st = (i-1)*code_param.symbols_per_frame+1;
+    en = st + code_param.symbols_per_frame-1;
+    llr(st:en) = sd_to_llr(r(st:en));    
+  end
+
+  % Outboard C decoder
+
+  if sdinput
+    f = fopen("sd.bin","wb"); fwrite(f, r, "double"); fclose(f);
+    system("../src/ldpc_dec sd.bin data_out.bin --sdinput"); 
+  else
+    f = fopen("llr.bin","wb"); fwrite(f, llr, "double"); fclose(f);
+    system("../src/ldpc_dec llr.bin data_out.bin"); 
+  end
+
+  f = fopen("data_out.bin","rb"); data_out = fread(f, "uint8")'; fclose(f);
+  
+  Nerrs = Nerrs2 = zeros(1,frames);;
+  for i=1:frames
+
+    % Check C decoder
+    
+    data_st = (i-1)*code_param.data_bits_per_frame+1;
+    data_en = data_st+code_param.data_bits_per_frame-1;
+    st = (i-1)*code_param.symbols_per_frame+1;
+    en = st+code_param.data_bits_per_frame-1;
+    data_out_c = data_out(st:en);
+    error_positions = xor(data_out_c, data(data_st:data_en));
+    Nerrs(i) = sum(error_positions);
+
+    % Octave decoder 
+
+    st = (i-1)*code_param.symbols_per_frame+1; en = st+code_param.symbols_per_frame-1;
+    [detected_data Niters] = ldpc_decode(r(st:en), code_param, max_iterations, decoder_type);
+    st = (i-1)*code_param.data_bits_per_frame+1; en = st+code_param.data_bits_per_frame-1;
+    data_out_octave = detected_data(1:code_param.data_bits_per_frame);
+    error_positions = xor(data_out_octave, data(st:en));
+    Nerrs2(i) = sum(error_positions);
+    %printf("%4d ", Niters);
+  end
+  printf("Errors per frame:\nC.....:");
+  for i=1:frames
+    printf("%4d ", Nerrs(i));
+  end
+  printf("\nOctave:");
+  for i=1:frames
+    printf("%4d ", Nerrs2(i));
+  end
+  printf("\n");
+
+end
+
+
+
 % Start simulation --------------------------------------------------------
 
 more off;
@@ -185,7 +266,7 @@ rand('state',1);
 
 % binary flags to run various demos, e.g. "15" to run 1 .. 8
 
-demos = 32;
+demos = 64;
 
 if bitand(demos,1)
   printf("simple_ut....\n");
@@ -230,7 +311,10 @@ if bitand(demos,16)
   printf("Nerrs = %d\n", Nerrs);
 end
 
-
 if bitand(demos,32)
   test_c_encoder;
 end
+
+if bitand(demos,64)
+  test_c_decoder;
+end