initial checkin of ldpc simulation files

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

diff --git a/codec2-dev/octave/ldpc.m b/codec2-dev/octave/ldpc.m
new file mode 100644 (file)
index 0000000..f30d688
--- /dev/null
+++ b/codec2-dev/octave/ldpc.m
@@ -0,0 +1,140 @@
+% ldpc.m\r
+% ldpc functions\r
+\r
+% LDPC demo;   Bill Cowley \r
+% Call the CML routines and simulate one set of SNRs.   See test_ldpc1.m\r
+%\r
+% sim_in the input parameter structure\r
+% sim_out contains BERs and other stats for each value of SNR\r
+% resfile is the result file\r
+%\r
+% 4/oct/2013:   edited to use the WiMax eIRA codes \r
+%               see 'help InitializeWiMaxLDPC' for parameter values \r
+\r
+1;\r
+\r
+function dummy\r
+endfunction\r
+\r
+function sim_out = ldpc_proc(sim_in, resfile)\r
+\r
+Eprob = sim_in.Eprob;  \r
+\r
+framesize = sim_in.framesize;\r
+rate      = sim_in.rate;\r
+modulation = sim_in.modulation;\r
+mod_order = sim_in.mod_order;\r
+mapping   = sim_in.mapping;\r
+\r
+\r
+Lim_Ferrs = sim_in.Lim_Ferrs;\r
+Ntrials   = sim_in.Ntrials;\r
+Esvec     = sim_in.Esvec;\r
+deb       = sim_in.deb;\r
+\r
+\r
+demod_type = 0;\r
+decoder_type = 0;\r
+max_iterations = 100;\r
+code_param.bits_per_symbol = log2(mod_order);\r
+bps = code_param.bits_per_symbol;\r
+\r
+[code_param.H_rows, code_param.H_cols, code_param.P_matrix] = InitializeWiMaxLDPC( rate, sim_in.framesize,  0 );\r
+\r
+code_param.data_bits_per_frame = length(code_param.H_cols) - length( code_param.P_matrix );\r
+\r
+code_param.S_matrix = CreateConstellation( modulation, mod_order, mapping );\r
+\r
+errfilename = '/home/david/codec2-dev/octave/mod_test_2000_poor_4dB.err';\r
+fin = fopen(errfilename, "rb");\r
+err = fread(fin,Inf, "short");\r
+length(err)\r
+Ntrials = floor(length(err)/framesize)\r
+\r
+for ne = 1:length(Esvec)\r
+    Es = Esvec(ne);\r
+    EsNo = 10^(Es/10);\r
+\r
+    \r
+    Terrs = 0;  Tbits =0;  Ferrs =0;\r
+    for nn = 1: Ntrials\r
+        \r
+        data = round( rand( 1, code_param.data_bits_per_frame ) );\r
+        \r
+        codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix );\r
+        code_param.code_bits_per_frame = length( codeword );\r
+        Nsymb = code_param.code_bits_per_frame/bps;\r
+       \r
+        % modulate\r
+        s = Modulate( codeword, code_param.S_matrix );\r
+        code_param.symbols_per_frame = length( s );\r
+     \r
+        s = Modulate(codeword, code_param.S_matrix );\r
+        code_param.symbols_per_frame = length( s );\r
+        \r
+      \r
+        variance = 1/(2*EsNo);\r
+        noise = sqrt(variance)*( randn(1,code_param.symbols_per_frame) + ...\r
+            j*randn(1,code_param.symbols_per_frame) );\r
+        a=ones(1, code_param.symbols_per_frame);  \r
+        r = a.*s + noise;\r
+\r
+        Nr = length(r); \r
+%        erasures = rand(1,Nr)<Eprob; \r
+%        r(erasures) = 0;            \r
+        \r
+        st = (nn-1)*Nr+1;\r
+        en = nn*Nr;\r
+        evec = err(st:en);\r
+        r(find(evec)) = 0;\r
+\r
+        symbol_likelihood = Demod2D( r, code_param.S_matrix, EsNo);\r
+         \r
+        % initialize the extrinsic decoder input\r
+        input_somap_c = zeros(1, code_param.code_bits_per_frame );\r
+        bit_likelihood = Somap( symbol_likelihood, demod_type, input_somap_c );\r
+        \r
+        input_decoder_c = bit_likelihood(1:code_param.code_bits_per_frame);\r
+        \r
+        x_hat= MpDecode( -input_decoder_c, code_param.H_rows, code_param.H_cols, ...\r
+            max_iterations, decoder_type, 1, 1);\r
+        detected_data = x_hat(max_iterations,:);\r
+        error_positions = xor( detected_data(1:code_param.data_bits_per_frame), data );\r
+        Nerrs = sum( error_positions);\r
+        \r
+        if Nerrs>0, fprintf(1,'x'),  else fprintf(1,'.'),  end\r
+        if (rem(nn, 50)==0),  fprintf(1,'\n'),  end\r
+        \r
+        \r
+        \r
+        if Nerrs>0,  Ferrs = Ferrs +1;  end\r
+        Terrs = Terrs + Nerrs;\r
+        Tbits = Tbits + code_param.data_bits_per_frame;\r
+        \r
+        if Ferrs > Lim_Ferrs, disp(['exit loop with #cw errors = ' ...\r
+                num2str(Ferrs)]);  break,  end\r
+    end\r
+    \r
+    TERvec(ne) = Terrs;\r
+    FERvec(ne) = Ferrs;\r
+    BERvec(ne) = Terrs/ Tbits;\r
+    Ebvec = Esvec - 10*log10(code_param.bits_per_symbol * rate);\r
+    \r
+    cparams= [code_param.data_bits_per_frame  code_param.symbols_per_frame ...\r
+        code_param.code_bits_per_frame];\r
+    \r
+    sim_out.BERvec = BERvec;\r
+    sim_out.Ebvec = Ebvec;\r
+    sim_out.FERvec = FERvec;\r
+    sim_out.TERvec  = TERvec;\r
+    sim_out.cpumins = cputime/60;\r
+    \r
+    if length(resfile)>0\r
+        save(resfile,  'sim_in',  'sim_out',  'cparams');\r
+        disp(['Saved results to ' resfile '  at Es =' num2str(Es) 'dB']);\r
+    end\r
+  end\r
+  \r
+endfunction\r
+\r
+\r

diff --git a/codec2-dev/octave/ldpcenc.m b/codec2-dev/octave/ldpcenc.m
new file mode 100644 (file)
index 0000000..233c9a3
--- /dev/null
+++ b/codec2-dev/octave/ldpcenc.m
@@ -0,0 +1,49 @@
+%test batch file   for the FSO-OOK simulations\r
+% \r
+% this version uses WiMax eIRA codes and includes erasures \r
+\r
+currentdir = pwd;\r
+thiscomp = computer;\r
+\r
+addpath '/home/david/tmp/cml/mat'    % assume the source files stored here\r
+cd /home/david/tmp/cml\r
+CmlStartup              % note that this is not in the cml path!\r
+disp('added cluster path and run CmlStartup')\r
+\r
+cd(currentdir)\r
+\r
+ldpc;\r
+\r
+%\r
+sim_in.Eprob = 0.1; \r
+disp([' test with erasure probability of ' num2str(sim_in.Eprob)]); \r
+\r
+sim_in.comment = 'test ldpc';\r
+sim_in.Esvec = 8:1/2:11; \r
+\r
+%sim_in.rate = 3/4;\r
+%sim_in.framesize = 16200\r
+\r
+% sim_in.rate = 0.5;\r
+% sim_in.framesize = 204;\r
+\r
+sim_in.rate = 3/4; \r
+sim_in.framesize = 576;  \r
+\r
+sim_in.mod_order = 4; \r
+sim_in.modulation = 'QPSK';\r
+sim_in.mapping = 'gray';\r
+\r
+sim_in.Lim_Ferrs= 30;\r
+sim_in.Ntrials =  1000;\r
+\r
+\r
+if exist('deb')~=1, deb = 0; end\r
+sim_in.deb =    deb;\r
+\r
+sim_out = ldpc_proc(sim_in, 'test.mat');\r
+\r
+figure(100); \r
+semilogy(sim_out.Ebvec,  sim_out.BERvec)\r
+xlabel('Eb/N0')\r
+ylabel('BER')\r