first pass at LDPC integration. Good results for AWGN channel, but HF fading not...
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 4 May 2017 05:36:43 +0000 (05:36 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 4 May 2017 05:36:43 +0000 (05:36 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@3122 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/ofdm_dev.m
codec2-dev/octave/ofdm_lib.m
codec2-dev/octave/ofdm_load_const.m

index 60377290f589651966ab36d5b3e26ab9f1d57439..b5ee2893d751aac006e4bffeff0328f277774851 100644 (file)
@@ -17,7 +17,8 @@ ofdm_lib;
         + plot error versus freq and timing offset
         + plot pro acquisition versus freq offset, timing and freq sep and together
         + plot total acquist prob at various SNRs ... maybe mean number of frames to sync?
-
+    [ ] replace genie EsNo est used for ldpc dec
+    [ ] clean up text 
 #}
 
 function [sim_out rx states] = run_sim(sim_in)
@@ -36,7 +37,7 @@ function [sim_out rx states] = run_sim(sim_in)
   timing_en = states.timing_en = sim_in.timing_en;
   states.foff_est_en = foff_est_en = sim_in.foff_est_en;
   states.phase_est_en = phase_est_en = sim_in.phase_est_en;
-
   if verbose
     printf("Rs:..........: %4.2f\n", Rs);
     printf("M:...........: %d\n", M);
@@ -46,7 +47,7 @@ function [sim_out rx states] = run_sim(sim_in)
     printf("Nrowsperframe: %d\n", Nrowsperframe);
     printf("Nsamperframe.: %d\n", Nsamperframe);
   end
-
+  
   % Important to define run time in seconds so HF model will evolve the same way
   % for different pilot insertion rates.  So lets work backwards from approx
   % seconds in run to get Nbits, the total number of payload data bits
@@ -73,6 +74,34 @@ function [sim_out rx states] = run_sim(sim_in)
     printf("Nrp....: %d (number of rows including pilots)\n", Nrp);
   end
 
+  % Optional LPDC code -----------------------------------------------
+
+  if isfield(sim_in, "ldpc_code") 
+    assert(bps == 2, "Only QPSK supported for LDPC so far.....");
+    HRA = sim_in.ldpc_code;
+    [aNr aNc] = size(HRA);
+    rate = states.rate = (aNc-aNr)/aNc;
+    assert(aNc == Nbitsperframe, "Dude: Num cols of LDPC HRA must == Nbitsperframe");
+    [H_rows, H_cols] = Mat2Hrows(HRA); 
+    code_param.H_rows = H_rows; 
+    code_param.H_cols = H_cols;
+    code_param.P_matrix = [];
+    code_param.data_bits_per_frame = length(code_param.H_cols) - length( code_param.P_matrix ); 
+    code_param.code_bits_per_frame = aNc;
+    assert(aNr == Nbitsperframe*rate);
+
+    ldpc_en = states.ldpc_en = 1;
+    modulation = states.ldpc_modulation = 'QPSK';
+    mapping = states.ldpc_mapping = 'gray';
+    demod_type = states.ldpc_demod_type = 0;
+    decoder_type = states.ldpc_decoder_type = 0;
+    max_iterations = states.ldpc_max_iterations = 100;
+
+    code_param.S_matrix = CreateConstellation( modulation, 4, mapping );
+
+    states.code_param = code_param;
+  end
+
   % set up HF model ---------------------------------------------------------------
 
   if hf_en
@@ -114,18 +143,26 @@ function [sim_out rx states] = run_sim(sim_in)
     rand('seed',1);
     randn('seed',1);
 
-    EbNo = bps * (10 .^ (EbNodB(nn)/10));
-    variance = 1/(M*EbNo/2);
+    EsNo = bps * (10 .^ (EbNodB(nn)/10));
+    variance = 1/(M*EsNo/2);
 
     Nsam = Nrp*(M+Ncp);
 
     % generate tx bits and run OFDM modulator
+    % note for reasons unknown LdpcEncode() returns garbage if we use > 0.5 rather than round()
 
-    tx_bits = rand(1,Nbits) > 0.5;
+    tx_data_bits = round(rand(1,Nbits*rate));
 
-    tx = [];
+    tx = []; tx_bits = [];
     for f=1:Nframes
-      tx = [tx ofdm_mod(states, tx_bits((f-1)*Nbitsperframe+1:f*Nbitsperframe))];
+      st = (f-1)*Nbitsperframe*rate+1; en = st + Nbitsperframe*rate - 1;
+      if ldpc_en
+        codeword = LdpcEncode(tx_data_bits(st:en), code_param.H_rows, code_param.P_matrix);
+      else
+        codeword = tx_data_bits(st:en);
+      end
+      tx = [tx ofdm_mod(states, codeword)];
+      tx_bits = [tx_bits codeword];
     end
 
     % add extra row of pilots at end, to allow one frame simulations,
@@ -179,7 +216,7 @@ function [sim_out rx states] = run_sim(sim_in)
     delta_t_log = []; 
     timing_est_log = [];
     foff_est_hz_log = [];
-    Nerrs_log = [];
+    Nerrs_log = []; Nerrs_coded_log = [];
     rx_bits = []; rx_np = [];
 
     % reset some states for each EbNo simulation point
@@ -225,13 +262,34 @@ function [sim_out rx states] = run_sim(sim_in)
     % calculate BER stats as a block, after pilots extracted
 
     errors = xor(tx_bits, rx_bits);
-    Nerrs = sum(errors);
+    Terrs = sum(errors);
+
+    Terrs_coded = 0; Tper_coded = 0;
     for f=1:Nframes
-      st = (f-1)*Nbitsperframe+1; en = st + Nbitsperframe-1;
+      st = (f-1)*Nbitsperframe+1; en = st + Nbitsperframe - 1;
       Nerrs_log(f) = sum(xor(tx_bits(st:en), rx_bits(st:en)));
+
+      % optional LDPC decode
+     
+      if ldpc_en
+        st = (f-1)*Nbitsperframe/bps + 1;
+        en = st + Nbitsperframe/bps - 1;
+        rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_np(st:en), EsNo);
+        st = (f-1)*Nbitsperframe*rate + 1;
+        en = st + Nbitsperframe*rate - 1;
+        Nerrs_coded = sum(xor(tx_data_bits(st:en), rx_codeword(1:Nbitsperframe*rate)));
+        Nerrs_coded_log(f) = Nerrs_coded;
+        Terrs_coded += Nerrs_coded;
+        if Nerrs_coded
+          Tper_coded++;
+        end
+      end
     end
 
-    printf("EbNodB: %3.2f BER: %5.4f Nbits: %d Nerrs: %d\n", EbNodB(nn), Nerrs/Nbits, Nbits, Nerrs);
+    printf("EbNodB: %3.2f BER: %5.4f Tbits: %d Terrs: %d\n", EbNodB(nn), Terrs/Nbits, Nbits, Terrs);
+    if ldpc_en
+      printf("       Coded BER: %5.4f Tbits: %d Terrs: %d PER: %5.4f\n", Terrs_coded/(Nbits*rate), Nbits*rate, Terrs_coded, Tper_coded/Nframes);
+    end
 
     if verbose
 
@@ -259,8 +317,17 @@ function [sim_out rx states] = run_sim(sim_in)
       title('Fine Freq');
 
       figure(5); clf;
-      plot(Nerrs_log);
-
+      if ldpc_en
+        subplot(211)
+        plot(Nerrs_log);
+        title("Uncoded Errors/frame");
+        subplot(212)
+        plot(Nerrs_coded_log);
+        title("Coded Errors/frame");
+     else
+        title("Errors/frame");
+        plot(Nerrs_log);
+      end
 
       figure(6)
       Tx = abs(fft(tx(1:Nsam).*hanning(Nsam)'));
@@ -293,7 +360,7 @@ function [sim_out rx states] = run_sim(sim_in)
 
     end
 
-    sim_out.ber(nn) = sum(Nerrs)/Nbits; 
+    sim_out.ber(nn) = Terrs/Nbits; 
     sim_out.pilot_overhead = 10*log10(Ns/(Ns-1));
   end
 endfunction
@@ -304,10 +371,10 @@ function run_single
   sim_in.Tcp = 0.002; 
   sim_in.Rs = 1/Ts; sim_in.bps = 2; sim_in.Nc = 16; sim_in.Ns = 8;
 
-  sim_in.Nsec = 5*(sim_in.Ns+1)/sim_in.Rs;  % one frame
-  %sim_in.Nsec = 30;
+  %sim_in.Nsec = 5*(sim_in.Ns+1)/sim_in.Rs;  % one frame
+  sim_in.Nsec = 30;
 
-  sim_in.EbNodB = 10;
+  sim_in.EbNodB = 100;
   sim_in.verbose = 1;
   sim_in.hf_en = 0;
   sim_in.foff_hz = 0;
@@ -317,6 +384,9 @@ function run_single
   sim_in.foff_est_en = 1;
   sim_in.phase_est_en = 1;
 
+  load HRA_112_112.txt
+  sim_in.ldpc_code = HRA_112_112;
+
   run_sim(sim_in);
 end
 
@@ -522,7 +592,9 @@ endfunction
 format;
 more off;
 
-%run_single
+init_cml('/home/david/Desktop/cml/');
+
+run_single
 %run_curves
 %acquisition_histograms
-acquisition_test
+%acquisition_test
index ebc3f28979908aed323ab1a1a914774801f341e4..3c1587c540e09cbe1e66a7e72f3ee7487072fe0a 100644 (file)
@@ -165,10 +165,36 @@ function states = ofdm_init(bps, Rs, Tcp, Ns, Nc)
 
   rate_fs_pilot_samples = states.pilots * W/states.M;
   states.rate_fs_pilot_samples = [rate_fs_pilot_samples(states.M-states.Ncp+1:states.M) rate_fs_pilot_samples];
+  
+  % LDPC code is optionally enabled
+
+  states.rate = 1.0;
+  states.ldpc_en = 0;
 
 endfunction
 
 
+% NOTE: You will need to set the CML path in the call to init_cml() below
+%       for you CML install.  See lpdc.m for instructions on how to install 
+%       CML library
+
+function init_cml(path_to_cml)
+  currentdir = pwd;
+  
+  if exist(path_to_cml, 'dir') == 7
+    cd(path_to_cml)
+    CmlStartup      
+    cd(currentdir); 
+  else
+    printf("\n---------------------------------------------------\n");
+    printf("Can't start CML in path: %s\n", path_to_cml);
+    printf("See CML path instructions at top of this script\n");
+    printf("-----------------------------------------------------\n\n");
+    assert(0);
+  end
+end
+
+
 % --------------------------------------
 % ofdm_mod - modulates one frame of bits
 % --------------------------------------
@@ -383,3 +409,22 @@ function [rx_bits states aphase_est_pilot_log rx_np] = ofdm_demod(states, rxbuf_
 endfunction
 
 
+function detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r, EsNo, fading)
+    if nargin == 6
+      fading = ones(1, length(r));
+    end
+
+    symbol_likelihood = Demod2D( r, code_param.S_matrix, EsNo, fading);
+         
+    % initialize the extrinsic decoder input
+
+    input_somap_c = zeros(1, code_param.code_bits_per_frame );
+    bit_likelihood = Somap( symbol_likelihood, demod_type, input_somap_c );
+        
+    input_decoder_c = bit_likelihood(1:code_param.code_bits_per_frame);
+        
+    x_hat= MpDecode( -input_decoder_c, code_param.H_rows, code_param.H_cols, ...
+                     max_iterations, decoder_type, 1, 1);
+    detected_data = x_hat(max_iterations,:);
+endfunction
+
index eb13e8e7c0d4e4a79fc01c31c04f878724fcfb0e..db66d6d98031fa8cd9d5f3fee2490361676b0386 100644 (file)
@@ -33,5 +33,14 @@ timing_en = states.timing_en;
 foff_est_en = states.foff_est_en;
 phase_est_en = states.phase_est_en;
 
+rate = states.rate;
+ldpc_en = states.ldpc_en;
+if ldpc_en
+  code_param = states.code_param;
+  max_iterations = states.ldpc_max_iterations;
+  demod_type = states.ldpc_demod_type;
+  decoder_type = states.ldpc_decoder_type;
+end
+
 verbose = states.verbose;