--- /dev/null
+% lpdc_fsk_lib.m\r
+% April 2015\r
+%\r
+% Library version of ldpc4.m written by vk5dsp. Application is high bit rate\r
+% balloon telemtry\r
+%\r
+% LDPC demo\r
+% Call the CML routines and simulate one set of SNRs. \r
+% This fucntion is an updated version of ldpc3() which uses less \r
+% of the CML functions \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
+\r
+1;\r
+\r
+function sim_out = ldpc5(sim_in, resfile, testmode, genie_Es, logging=0);\r
+\r
+ if nargin<4, testmode = 0; end\r
+ estEsN0 = 0; \r
+\r
+ HRA = sim_in.HRA;\r
+ framesize = sim_in.framesize;\r
+ rate = sim_in.rate;\r
+ mod_order = sim_in.mod_order;\r
+\r
+ Lim_Ferrs = sim_in.Lim_Ferrs;\r
+ Ntrials = sim_in.Ntrials;\r
+ Esvec = sim_in.Esvec;\r
+\r
+ demod_type = 0;\r
+ decoder_type = 0;\r
+ max_iterations = 100;\r
+ code_param = ldpc_init(HRA, mod_order);\r
+ bps = code_param.bits_per_symbol;\r
+\r
+\r
+ if (logging) \r
+ fod = fopen('decode.log', 'w'); \r
+ fwrite(fod, 'Es estEs Its secs \n'); \r
+ end \r
+\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
+ codeword = ldpc_encode(code_param, data);\r
+\r
+ code_param.code_bits_per_frame = length( codeword );\r
+ Nsymb = code_param.code_bits_per_frame/bps;\r
+\r
+ if testmode==1\r
+ f1 = fopen("dat_in2064.txt", "w");\r
+ for k=1:length(data); fprintf(f1, "%u\n", data(k)); end\r
+ fclose(f1); \r
+ system("./ra_enc"); \r
+\r
+ load("dat_op2064.txt"); \r
+ pbits = codeword(length(data)+1:end); % print these to compare with C code \r
+ dat_op2064(1:16)', pbits(1:16) \r
+ differences_in_parity = sum(abs(pbits - dat_op2064'))\r
+ pause; \r
+ end\r
+\r
+\r
+ % modulate\r
+ % s = Modulate( codeword, code_param.S_matrix );\r
+ s= 1 - 2 * codeword; \r
+ code_param.symbols_per_frame = length( s );\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
+ r = s + noise;\r
+ Nr = length(r); \r
+\r
+ [detected_data Niters] = ldpc_decode(r, code_param, max_iterations, decoder_type);\r
+\r
+ error_positions = xor( detected_data(1:code_param.data_bits_per_frame), data );\r
+ Nerrs = sum( error_positions);\r
+\r
+ t = clock; t = fix(t(5)*60+t(6)); \r
+ if (logging) \r
+ fprintf(fod, ' %3d %4d\n', Niters, t); \r
+ end \r
+\r
+ if Nerrs>0, fprintf(1,'x'), else fprintf(1,'.'), end\r
+ if (rem(nn, 50)==0), fprintf(1,'\n'), end\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 nargin > 2\r
+ save(resfile, 'sim_in', 'sim_out', 'cparams');\r
+ disp(['Saved results to ' resfile ' at Es =' num2str(Es) 'dB']);\r
+ end\r
+ end\r
+end\r
+\r
+\r
+function code_param = ldpc_init(HRA, mod_order)\r
+ code_param.bits_per_symbol = log2(mod_order);\r
+ [H_rows, H_cols] = Mat2Hrows(HRA); \r
+ code_param.H_rows = H_rows; \r
+ code_param.H_cols = H_cols;\r
+ code_param.P_matrix = [];\r
+ code_param.data_bits_per_frame = length(code_param.H_cols) - length( code_param.P_matrix ); \r
+ code_param.symbols_per_frame = length(HRA);\r
+end\r
+\r
+\r
+function codeword = ldpc_encode(code_param, data)\r
+ codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix );\r
+endfunction\r
+\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 / 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
+ [x_hat, PCcnt] = MpDecode( input_decoder_c, 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
+ 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
+ end\r
+end\r
+\r
+\r
+% optionally create a C include file for use in mpdecode.c C cmd line LDPC decoder\r
+\r
+function write_code_to_C_include_file(code_param, max_iterations, decoder_type, input_decoder_c, x_hat)\r
+ \r
+ f = fopen(code_param.c_include_file, "wt");\r
+\r
+ fprintf(f, "/*\n FILE....: %s\n\n Static arrays for TRC HF modem, generated", code_param.c_include_file);\r
+ fprintf(f, "\n by test_ldpc_fsk.m:simple_ut().\n\n*/\n\n");\r
+ \r
+ fprintf(f,"#define NUMBERPARITYBITS %d\n", rows(code_param.H_rows));\r
+ fprintf(f,"#define MAX_ROW_WEIGHT %d\n", columns(code_param.H_rows));\r
+ fprintf(f,"#define CODELENGTH %d\n", code_param.symbols_per_frame);\r
+ fprintf(f,"#define NUMBERROWSHCOLS %d\n", rows(code_param.H_cols));\r
+ fprintf(f,"#define MAX_COL_WEIGHT %d\n", columns(code_param.H_cols));\r
+ fprintf(f,"#define DEC_TYPE %d\n", decoder_type);\r
+ fprintf(f,"#define MAX_ITER %d\n", max_iterations);\r
+\r
+ fprintf(f,"\ndouble H_rows[] = {\n");\r
+ \r
+ % clock out 2D array to linear C array in row order ....\r
+\r
+ [r c] = size(code_param.H_rows);\r
+ for j=1:c\r
+ for i=1:r\r
+ fprintf(f, "%f", code_param.H_rows(i,j));\r
+ if (i == r) && (j ==c)\r
+ fprintf(f,"\n};\n");\r
+ else\r
+ fprintf(f,", ");\r
+ end\r
+ end\r
+ end\r
+\r
+ fprintf(f,"\ndouble H_cols[] = {\n");\r
+ [r c] = size(code_param.H_cols);\r
+ for j=1:c\r
+ for i=1:r\r
+ fprintf(f, "%f", code_param.H_cols(i,j));\r
+ if (i == r) && (j == c)\r
+ fprintf(f,"\n};\n");\r
+ else\r
+ fprintf(f,", ");\r
+ end\r
+ end\r
+ end\r
+\r
+ fprintf(f,"\ndouble input[] = {\n");\r
+ for i=1:length(input_decoder_c)\r
+ fprintf(f, "%f", input_decoder_c(i));\r
+ if i == length(input_decoder_c)\r
+ fprintf(f,"\n};\n");\r
+ else\r
+ fprintf(f,", "); \r
+ end\r
+ end\r
+\r
+ fprintf(f,"\nint output[] = {\n");\r
+ [r c] = size(x_hat);\r
+ for j=1:c\r
+ for i=1:r\r
+ fprintf(f, "%d", x_hat(i,j));\r
+ if (i == r) && (j == c)\r
+ fprintf(f,"\n};\n");\r
+ else\r
+ fprintf(f,", ");\r
+ end\r
+ end\r
+ end\r
+\r
+ fclose(f);\r
+end\r
--- /dev/null
+% test_ldpc_fsk_lib\r
+% David Rowe 16 April 2016\r
+%\r
+% Some tests for ldpc_fsk_lib\r
+\r
+1;\r
+\r
+% encodes and decodes one frame, also writes codeword.bin for testing\r
+% decode_from_file() below, and can optionally generate include file for\r
+% C version of decoder.\r
+\r
+function data = simple_ut(c_include_file)\r
+ load('H2064_516_sparse.mat');\r
+ HRA = full(HRA); \r
+ max_iterations = 100;\r
+ decoder_type = 0;\r
+ EsNodB = 4;\r
+ mod_order = 2;\r
+\r
+ code_param = ldpc_init(HRA, mod_order);\r
+ data = round( rand( 1, code_param.data_bits_per_frame ) );\r
+ codeword = ldpc_encode(code_param, data);\r
+ f = fopen("codeword.bin","wt"); fwrite(f, codeword, "uint8"); fclose(f);\r
+ s = 1 - 2 * codeword; \r
+ code_param.symbols_per_frame = length( s );\r
+\r
+ EsNo = 10^(EsNodB/10);\r
+ variance = 1/(2*EsNo);\r
+ noise = sqrt(variance)* randn(1,code_param.symbols_per_frame); \r
+ rx = s + noise;\r
+ \r
+ if nargin == 1\r
+ code_param.c_include_file = c_include_file;\r
+ end\r
+ [detected_data Niters] = ldpc_decode(rx, code_param, max_iterations, decoder_type);\r
+ \r
+ error_positions = xor(detected_data(1:code_param.data_bits_per_frame), data);\r
+ Nerrs = sum(error_positions);\r
+\r
+ printf(" Nerrs = %d\n", Nerrs);\r
+\r
+end\r
+\r
+\r
+% This version decodes from a file of bits\r
+\r
+function detected_data = decode_from_file(filename)\r
+ max_iterations = 100;\r
+ decoder_type = 0;\r
+ load('H2064_516_sparse.mat');\r
+ HRA = full(HRA); \r
+ mod_order = 2;\r
+\r
+ f = fopen(filename,"rb"); codeword = fread(f, "uint8")'; fclose(f);\r
+ r = 1 - 2 * codeword; \r
+ code_param = ldpc_init(HRA, mod_order);\r
+ [detected_data Niters] = ldpc_decode(r, code_param, max_iterations, decoder_type);\r
+end\r
+\r
+\r
+% plots a BER cuvre for the decoder. Takes a while to run, uses parallel cores\r
+\r
+function plot_curve\r
+ num_cores = 4; % set this to the number of cores you have\r
+\r
+ load('H2064_516_sparse.mat');\r
+ HRA = full(HRA); \r
+ [Nr Nc] = size(HRA); \r
+ sim_in.rate = (Nc-Nr)/Nc;\r
+\r
+ sim_in.HRA = HRA;\r
+ sim_in.mod_order = 2;\r
+ sim_in.framesize = Nc;\r
+ sim_in.mod_order = 2; \r
+ sim_in.Lim_Ferrs = 100;\r
+\r
+ % note we increase number of trials as BER goes down\r
+\r
+ Esvec = [ 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 ]; \r
+ Ntrials = [ 1E4 1E4 1E4 1E4 1E5 1E5 1E5 1E5 1E5 ];\r
+ num_runs = length(Esvec)\r
+\r
+ sim_in_vec(1:num_runs) = sim_in;\r
+ for i = 1:num_runs\r
+ sim_in_vec(i).Esvec = Esvec(i);\r
+ sim_in_vec(i).Ntrials = Ntrials(i);\r
+ end\r
+\r
+ %sim_out = ldpc5(sim_in_vec(1));\r
+ tstart = time();\r
+ sim_out = pararrayfun(num_cores, @ldpc5, sim_in_vec);\r
+ tend = time();\r
+\r
+ total_bits = sum(Ntrials)*sim_in.framesize;\r
+ total_secs = tend - tstart;\r
+ printf("%d bits in %4.1f secs, or %5f bits/s\n", total_bits, total_secs, total_bits/total_secs);\r
+\r
+ for i=1:num_runs\r
+ Ebvec(i) = sim_out(i).Ebvec;\r
+ BERvec(i) = sim_out(i).BERvec;\r
+ end\r
+ semilogy(Ebvec, BERvec, '+-')\r
+ xlabel('Eb/N0')\r
+ ylabel('BER')\r
+ title(['H2064 516 sparse.mat' ' ' date])\r
+\r
+end\r
+\r
+\r
+% Start simulation --------------------------------------------------------\r
+\r
+more off;\r
+currentdir = pwd;\r
+thiscomp = computer;\r
+\r
+if strcmpi(thiscomp, 'MACI64')==1\r
+ if exist('CMLSimulate')==0\r
+ cd '/Users/bill/Current/Projects/DLR_FSO/Visit2013_FSO_GEO/cml'\r
+ addpath '../' % assume the source files stored here\r
+ CmlStartup % note that this is not in the cml path!\r
+ disp('added MACI64 path and run CmlStartup')\r
+ end\r
+end\r
+\r
+if strfind(thiscomp, 'pc-linux-gnu')==8 \r
+ if exist('LdpcEncode')==0, \r
+ cd '~/tmp/cml'\r
+ CmlStartup \r
+ disp('CmlStartup has been run')\r
+ % rmpath '/home/bill/cml/mexhelp' % why is this needed? \r
+ % maybe different path order in octave cf matlab ? \r
+ end\r
+end\r
+\r
+cd(currentdir)\r
+\r
+ldpc_fsk_lib;\r
+randn('state',1);\r
+rand('state',1);\r
+\r
+% binary flags to run various demos, e.g. "15" runs them all\r
+\r
+demos = 2;\r
+\r
+if bitand(demos,1)\r
+ printf("simple_ut....\n");\r
+ data = simple_ut;\r
+end\r
+\r
+if bitand(demos,2)\r
+ printf("generate C header file....\n");\r
+ data = simple_ut("../src/ldpc_code.h");\r
+end\r
+\r
+if bitand(demos,4)\r
+ printf("decode_from_file ......\n");\r
+ detected_data = decode_from_file("codeword.bin");\r
+ error_positions = xor( detected_data(1:length(data)), data );\r
+ Nerrs = sum(error_positions);\r
+ printf(" Nerrs = %d\n", Nerrs);\r
+end\r
+\r
+if bitand(demos,8)\r
+ printf("plot a curve....\n");\r
+ plot_curve;\r
+end\r
+\r
projects / freetel-svn-tracking.git / commitdiff