diff options
Diffstat (limited to 'octave/ldpc.m')
| -rw-r--r-- | octave/ldpc.m | 194 |
1 files changed, 194 insertions, 0 deletions
diff --git a/octave/ldpc.m b/octave/ldpc.m new file mode 100644 index 0000000..e9d6ede --- /dev/null +++ b/octave/ldpc.m @@ -0,0 +1,194 @@ +% ldpc.m +% +#{ + David Rowe 2013 + Octave functions for the CML LDPC library. + + To install and compile CML support: + + $ sudo apt-get install liboctave-dev + $ git clone [email protected]:drowe67/cml.git + $ cd cml + $ make + + If you have configured codec2 with cmake -DUNITTEST=1, then you will + already have CML (e.g. under build_linux/cml), as it is used to run unit tests. + + To use CML when running Octave simulations from the Octave CLI, set an + environment variable for CML_PATH in your shell or in your + codec2/octave/.octaverc file: + + setenv("CML_PATH",sprintf("%s/codec2/build_linux/cml",getenv("HOME"))) +#} + +1; + +function init_cml() + currentdir = pwd; + + path_to_cml = getenv("CML_PATH"); + + 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 (ldpc.m)\n"); + printf("-----------------------------------------------------\n\n"); + assert(0); + end +end + +% init using built in WiMax or DVSB2 code + +function code_param = ldpc_init_builtin(code, rate, framesize, modulation, mod_order, mapping, constellation) + if strcmp(code,'wimax') + [code_param.H_rows, code_param.H_cols, code_param.P_matrix] = InitializeWiMaxLDPC( rate, framesize, 0 ); + end + if strcmp(code,'dvbs2') + [code_param.H_rows, code_param.H_cols, code_param.P_matrix] = InitializeDVBS2( rate, framesize); + end + if nargin == 7 + code_param.S_matrix = constellation; + else + if length(mapping) == 0 + code_param.S_matrix = CreateConstellation( modulation, mod_order); + else + code_param.S_matrix = CreateConstellation( modulation, mod_order, mapping ); + end + end + code_param.bits_per_symbol = log2(mod_order); + + code_param.ldpc_data_bits_per_frame = length(code_param.H_cols) - length(code_param.P_matrix); + code_param.ldpc_parity_bits_per_frame = framesize - code_param.ldpc_data_bits_per_frame; + code_param.ldpc_coded_bits_per_frame = framesize; + + code_param.data_bits_per_frame = code_param.ldpc_data_bits_per_frame; + code_param.coded_bits_per_frame = code_param.ldpc_coded_bits_per_frame; + code_param.coded_syms_per_frame = code_param.coded_bits_per_frame/code_param.bits_per_symbol; +endfunction + + +% init using user supplied code + +function [code_param framesize rate] = ldpc_init_user(HRA, modulation, mod_order, mapping, constellation) + [Nr Nc] = size(HRA); + rate = (Nc-Nr)/Nc; + framesize = Nc; + [H_rows, H_cols] = Mat2Hrows(HRA); + code_param.H_rows = H_rows; + code_param.H_cols = H_cols; + code_param.P_matrix = []; + if nargin == 5 + code_param.S_matrix = constellation; + else + if length(mapping) == 0 + code_param.S_matrix = CreateConstellation( modulation, mod_order); + else + code_param.S_matrix = CreateConstellation( modulation, mod_order, mapping ); + end + end + code_param.bits_per_symbol = log2(mod_order); + + code_param.ldpc_data_bits_per_frame = length(code_param.H_cols) - length(code_param.P_matrix); + code_param.ldpc_parity_bits_per_frame = framesize - code_param.ldpc_data_bits_per_frame; + code_param.ldpc_coded_bits_per_frame = framesize; + + % these variables support 1's stuffing (not using all data bits to lower code rate) + code_param.data_bits_per_frame = code_param.ldpc_data_bits_per_frame; + code_param.coded_bits_per_frame = code_param.ldpc_coded_bits_per_frame; + code_param.coded_syms_per_frame = code_param.coded_bits_per_frame/code_param.bits_per_symbol; +endfunction + + +function [codeword s] = ldpc_enc(data, code_param) + if code_param.data_bits_per_frame != code_param.ldpc_data_bits_per_frame + % optionally lower the code rate by "1's stuffing" - setting Nunused data bits to 1 + Nunused = code_param.ldpc_data_bits_per_frame - code_param.data_bits_per_frame; + codeword = LdpcEncode([data ones(1,Nunused)], code_param.H_rows, code_param.P_matrix); + % remove unused data bits from codeword, as they are known to the receiver and don't need to be transmitted + codeword = [ codeword(1:code_param.data_bits_per_frame) codeword(code_param.ldpc_data_bits_per_frame+1:end) ]; + else + codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix ); + end + s = Modulate( codeword, code_param.S_matrix ); +endfunction + + +function [detected_data paritychecks] = ldpc_dec(code_param, max_iterations, ... + demod_type, decoder_type, r, ... + EsNo, fading) + % handle "1's stuffing" case where we don't use all data bits + Nunused = code_param.ldpc_data_bits_per_frame - code_param.data_bits_per_frame; + + symbol_likelihood = Demod2D( r, code_param.S_matrix, EsNo, fading); + + % initialize the extrinsic decoder input + + input_somap_c = zeros(1, code_param.ldpc_coded_bits_per_frame - Nunused); + bit_likelihood = Somap( symbol_likelihood, demod_type, input_somap_c ); + + input_decoder_c = bit_likelihood(1:(code_param.ldpc_coded_bits_per_frame-Nunused)); + + % insert "very likely" LLRs for unused data bits (in 1's stuffing case) + input_decoder_c = [input_decoder_c(1:code_param.data_bits_per_frame) ... + 100*ones(1,Nunused) ... + input_decoder_c(code_param.data_bits_per_frame+1:end)]; + + [x_hat paritychecks] = MpDecode( -input_decoder_c, code_param.H_rows, code_param.H_cols, ... + max_iterations, decoder_type, 1, 1); + [mx mx_ind] = max(paritychecks); + detected_data = x_hat(mx_ind,:); +endfunction + + +% Packs a binary array into an array of 8 bit bytes, MSB first + +function packed = packmsb(unpacked) + packed = zeros(1,floor(length(unpacked)+7)/8); + bit = 7; byte = 1; + for i=1:length(unpacked) + packed(byte) = bitor(packed(byte), bitshift(unpacked(i),bit)); + bit--; + if (bit < 0) + bit = 7; + byte++; + end + end +endfunction + + +% unpacks an array of 8 bit bytes into a binary array of unpacked bits, MSB first + +function unpacked = unpackmsb(packed) + bit = 7; byte = 1; + for i=1:length(packed)*8 + unpacked(i) = bitand(bitshift(packed(byte), -bit), 1); + bit--; + if (bit < 0) + bit = 7; + byte++; + end + end +endfunction + + +% symbol interleaver that acts on bits 2 at a time + +function y = interleave_bits(interleaver, x) + y = zeros(1,length(x)); + for i = 1:length(interleaver) + dst = interleaver(i); + y(2*(dst-1)+1:2*dst) = x(2*(i-1)+1:2*(i)); + end +endfunction + +% symbol de-interleaver + +function x = deinterleave_symbols(interleaver, y) + for i = 1:length(interleaver) + x(i) = y(interleaver(i)); + end +endfunction |