3 step ldpc demo withgd results

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

diff --git a/codec2-dev/octave/ldpcut.m b/codec2-dev/octave/ldpcut.m
index 39843d8d7d373e7845121f46555af5a2a6cc2917..f6d787ef7da80f6fd1509910ee6c36a00a205595 100644 (file)
--- a/codec2-dev/octave/ldpcut.m
+++ b/codec2-dev/octave/ldpcut.m
@@ -17,58 +17,141 @@ cd(currentdir)
 \r
 ldpc;\r
 \r
-% Start simulation\r
+function sim_out = run_simulation(sim_in)\r
 \r
-rate = 3/4; \r
-framesize = 576;  \r
+  rate = sim_in.rate; \r
+  framesize = sim_in.framesize;  \r
+  Ntrials = sim_in.Ntrials;\r
+  EsNodBvec = sim_in.EsNodBvec;\r
+  verbose = sim_in.verbose;\r
+\r
+  % Start simulation\r
+\r
+  mod_order = 4; \r
+  modulation = 'QPSK';\r
+  mapping = 'gray';\r
+\r
+  demod_type = 0;\r
+  decoder_type = 0;\r
+  max_iterations = 100;\r
+\r
+  code_param = ldpc_init(rate, framesize, modulation, mod_order, mapping);\r
+\r
+  for ne = 1:length(EsNodBvec)\r
+    EsNodB = EsNodBvec(ne);\r
+    EsNo = 10^(EsNodB/10);\r
+    variance = 1/EsNo;\r
+\r
+    Tbits = Terrs = Ferrs = 0;\r
+    \r
+    tx_bits = [];\r
+    tx_symbols = []; \r
+    rx_symbols = [];\r
+\r
+    % Encode a bunch of frames\r
+\r
+    for nn = 1: Ntrials        \r
+      atx_bits = round(rand( 1, code_param.data_bits_per_frame));\r
+      tx_bits = [tx_bits atx_bits];\r
+      [tx_codeword atx_symbols] = ldpc_enc(atx_bits, code_param);\r
+      tx_symbols = [tx_symbols atx_symbols];\r
+    end\r
+\r
+    % Add AWGN noise, 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im\r
+\r
+    noise = sqrt(variance*0.5)*(randn(1,length(tx_symbols)) + j*randn(1,length(tx_symbols)));\r
+    rx_symbols = tx_symbols + noise;\r
+\r
+    % Decode a bunch of frames\r
+\r
+    for nn = 1: Ntrials        \r
+      st = (nn-1)*code_param.symbols_per_frame + 1;\r
+      en = (nn)*code_param.symbols_per_frame;\r
+      arx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symbols(st:en), EsNo, ones(1,code_param.symbols_per_frame));\r
+      st = (nn-1)*code_param.data_bits_per_frame + 1;\r
+      en = (nn)*code_param.data_bits_per_frame;\r
+      error_positions = xor(arx_codeword(1:code_param.data_bits_per_frame), tx_bits(st:en));\r
+      Nerrs = sum( error_positions);\r
+        \r
+      if verbose == 2\r
+        % print "." if frame decoded without errors, 'x' if we can't decode\r
+\r
+        if Nerrs > 0, printf('x'),  else printf('.'),  end\r
+        if (rem(nn, 50)==0),  printf('\n'),  end    \r
+      end\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
+    if verbose\r
+      printf("\nEsNodB: %3.2f BER: %f Tbits: %d Terrs: %d Ferrs: %d\n", EsNodB, Terrs/Tbits, Tbits, Terrs,  Ferrs)\r
+    end\r
+\r
+    sim_out.Tbits(ne) = Tbits;\r
+    sim_out.Terrs(ne) = Terrs;\r
+    sim_out.Ferrs(ne) = Ferrs;\r
+    sim_out.BER(ne)   = Terrs/Tbits;\r
+    sim_out.FER(ne)   = Ferrs/Ntrials;\r
+  end\r
+\r
+endfunction\r
+\r
+% START SIMULATIONS ---------------------------------------------------------------\r
+\r
+more off;\r
+\r
+% 1/ Simplest possible one frame simulation ---------------------------------------\r
+\r
+printf("Test 1\n------\n");\r
 \r
 mod_order = 4; \r
 modulation = 'QPSK';\r
 mapping = 'gray';\r
-\r
+framesize = 576;         % CML library has a bunch of different framesizes available\r
+rate = 1/2;\r
 demod_type = 0;\r
 decoder_type = 0;\r
 max_iterations = 100;\r
+EsNo = 10;               % decoder needs an estimated channel EsNo (linear ratio, not dB)\r
 \r
 code_param = ldpc_init(rate, framesize, modulation, mod_order, mapping);\r
+tx_bits = round(rand(1, code_param.data_bits_per_frame));\r
+[tx_codeword, qpsk_symbols] = ldpc_enc(tx_bits, code_param);\r
+rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, qpsk_symbols, EsNo);\r
 \r
-Ntrials = 84;\r
-EsNo = 10;\r
+errors_positions = xor(tx_bits, rx_codeword(1:framesize*rate));\r
+Nerr = sum(errors_positions);\r
+printf("Nerr: %d\n\n", Nerr);\r
 \r
-Tbits = Terrs = Ferrs = 0;\r
-    \r
-data = [];\r
-r = []; \r
-\r
-% Encode a bunch of frames\r
-\r
-for nn = 1: Ntrials        \r
-    d = round( rand( 1, code_param.data_bits_per_frame ) );\r
-    data = [data d];\r
-    [codeword, s] = ldpc_enc(d, code_param);\r
-    code_param.code_bits_per_frame = length(codeword);\r
-    code_param.symbols_per_frame = length(s);\r
-    r = [r s];\r
-end\r
-\r
-% Decode a bunch of frames\r
-\r
-for nn = 1: Ntrials        \r
-    st = (nn-1)*code_param.symbols_per_frame + 1;\r
-    en = (nn)*code_param.symbols_per_frame;\r
-    detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r(st:en), EsNo, ones(1,code_param.symbols_per_frame));\r
-    st = (nn-1)*code_param.data_bits_per_frame + 1;\r
-    en = (nn)*code_param.data_bits_per_frame;\r
-    error_positions = xor( detected_data(1:code_param.data_bits_per_frame), data(st:en) );\r
-    Nerrs = sum( error_positions);\r
-        \r
-    % print "." if frame decoded without errors, 'x' if we can't decode\r
+% 2/ Run a bunch of trials at just one EsNo point --------------------------------------\r
+\r
+printf("Test 2\n------\n");\r
+\r
+sim_in.rate = 0.5; \r
+sim_in.framesize = 2*576;  \r
+sim_in.verbose = 2;\r
+sim_in.Ntrials = 100;\r
+sim_in.EsNodBvec = 5;\r
+run_simulation(sim_in);\r
+\r
+% 3/ Lets draw a Eb/No versus BER curve -------------------------------------------------\r
+\r
+printf("\n\nTest 3\n------\n");\r
+\r
+sim_in.EsNodBvec = -2:10;\r
+sim_out = run_simulation(sim_in);\r
 \r
-    if Nerrs>0, fprintf(1,'x'),  else printf('.'),  end\r
-    if (rem(nn, 50)==0),  printf('\n'),  end    \r
-    if Nerrs>0,  Ferrs = Ferrs +1;  end\r
-    Terrs = Terrs + Nerrs;\r
-    Tbits = Tbits + code_param.data_bits_per_frame;        \r
-end\r
-printf("\nTbits: %d Terrs: %d Ferrs: %d\n", Tbits, Terrs,  Ferrs)\r
+EbNodB = sim_in.EsNodBvec - 10*log10(2);  % QPSK has two bits/symbol\r
+uncoded_BER_theory = 0.5*erfc(sqrt(10.^(EbNodB/10)));\r
 \r
+figure(1)\r
+clf\r
+semilogy(EbNodB, uncoded_BER_theory,'r;uncoded QPSK theory;')\r
+hold on;\r
+semilogy(EbNodB-10*log10(sim_in.rate), sim_out.BER,'g;LDPC coded QPSK simulation;');\r
+hold off;\r
+grid('minor')\r
+xlabel('Eb/No (dB)')\r
+ylabel('BER')\r