LDPC codes with belief propagation decoding
0
fork

Configure Feed

Select the types of activity you want to include in your feed.

Remove self-test interop (LDPC, Turbo); extend SCC dariol83

LDPC and Turbo interop tests reimplemented the CCSDS encoder in
Python ("we implement the full CCSDS H-matrix construction",
"we implement the full CCSDS turbo encoder here"). Not interop.

SCC dariol83 interop extended with Reed-Solomon traces using
dariol83's RS codec.

-378
-17
test/interop/commpy/dune
··· 1 - (test 2 - (name test) 3 - (libraries ldpc csvt alcotest) 4 - (deps 5 - (source_tree traces) 6 - (source_tree scripts))) 7 - 8 - ; Regenerate traces against commpy: dune build @regen-traces 9 - 10 - (rule 11 - (alias regen-traces) 12 - (deps 13 - (source_tree scripts)) 14 - (action 15 - (chdir 16 - scripts 17 - (run bash ./generate.sh))))
-258
test/interop/commpy/scripts/generate.py
··· 1 - #!/usr/bin/env python3 2 - """Generate CCSDS LDPC interop traces for ocaml-ldpc. 3 - 4 - Oracle: standalone CCSDS AR4JA encoder 5 - Regenerate: dune build @regen-traces 6 - 7 - Implements the CCSDS 131.0-B-4 Section 7 AR4JA rate 1/2 LDPC encoder: 8 - - k=1024 information bits, n=2048 transmitted bits 9 - - n_full=2560 (512 punctured bits) 10 - - Circulant-based parity-check matrix with two-step expansion 11 - - Systematic encoding via Gaussian elimination over GF(2) 12 - 13 - commpy provides a generic LDPC BP decoder/encoder but does NOT include 14 - CCSDS AR4JA code constructions, so we implement the full CCSDS H-matrix 15 - construction and systematic encoder here, matching the OCaml implementation. 16 - """ 17 - import os 18 - import sys 19 - 20 - import numpy as np 21 - 22 - 23 - # --- Constants --- 24 - 25 - M_SUB = 512 # submatrix size 26 - CIRC = M_SUB // 4 # circulant size = 128 27 - NUM_ROWS = 3 # prototype check rows 28 - NUM_COLS = 5 # prototype variable columns (only first 5 non-zero) 29 - H_M = NUM_ROWS * M_SUB # 1536 check equations 30 - H_N = NUM_COLS * M_SUB # 2560 full codeword length 31 - K = 1024 # information bits 32 - N = 2048 # transmitted bits 33 - PUNCTURED = 512 # punctured bits 34 - 35 - HZ = 0 # zero block 36 - HI = 1 # identity (shifted) 37 - HP = 2 # permutation block 38 - 39 - # Prototype layers (3 rows x 11 columns, but only first 5 cols matter) 40 - # Layer 0 41 - PROTO_0 = [ 42 - [(HZ,0), (HZ,0), (HI,0), (HZ,0), (HI,0)], 43 - [(HI,0), (HI,0), (HZ,0), (HI,0), (HP,1)], 44 - [(HI,0), (HP,4), (HZ,0), (HP,6), (HI,0)], 45 - ] 46 - 47 - # Layer 1 48 - PROTO_1 = [ 49 - [(0,0), (0,0), (0,0), (0,0), (HP,0)], 50 - [(0,0), (0,0), (0,0), (0,0), (HP,2)], 51 - [(0,0), (HP,5), (0,0), (HP,7), (0,0)], 52 - ] 53 - 54 - # Layer 2 55 - PROTO_2 = [ 56 - [(0,0), (0,0), (0,0), (0,0), (0,0)], 57 - [(0,0), (0,0), (0,0), (0,0), (HP,3)], 58 - [(0,0), (0,0), (0,0), (0,0), (0,0)], 59 - ] 60 - 61 - # THETA_K: row permutation index 62 - THETA_K = [ 63 - 3, 0, 1, 2, 2, 3, 0, 1, 0, 1, 2, 0, 2, 3, 0, 1, 64 - 2, 0, 1, 2, 0, 1, 2, 1, 2, 3, 65 - ] 66 - 67 - # PHI_J_K for M=128 (circulant_size = 32, values mod 32) 68 - PHI_J_K = [ 69 - [1,22,0,26,0,10,5,18,3,22,3,8,25,25,2,27,7,7,15,10,4,19,7,9,26,17], 70 - [0,27,30,28,7,1,8,20,26,24,4,12,23,15,15,22,31,3,29,21,2,5,11,26,9,17], 71 - [0,12,30,18,10,16,13,9,7,15,16,18,4,23,5,3,29,11,4,8,2,11,11,3,15,13], 72 - [0,13,19,14,15,20,17,4,4,11,17,20,8,22,19,15,5,21,17,9,20,18,31,13,2,18], 73 - ] 74 - 75 - 76 - def build_ar4ja_h(): 77 - """Construct the CCSDS AR4JA parity-check matrix H for rate 1/2, k=1024. 78 - 79 - Returns a sparse representation: list of (row, col) pairs. 80 - """ 81 - edges = [[] for _ in range(H_M)] 82 - mod_m = M_SUB - 1 83 - mod_circ = CIRC - 1 84 - log_circ = 7 # log2(128) 85 - 86 - def process_layer(proto): 87 - for row_idx in range(NUM_ROWS): 88 - for col_idx in range(NUM_COLS): 89 - btype, bval = proto[row_idx][col_idx] 90 - if btype == HI: 91 - for c in range(M_SUB): 92 - chk = row_idx * M_SUB + c 93 - var = col_idx * M_SUB + ((c + bval) & mod_m) 94 - edges[chk].append(var) 95 - elif btype == HP: 96 - for c in range(M_SUB): 97 - j = c >> log_circ # j = c // 128 98 - pi = (((THETA_K[bval] + j) % 4) << log_circ) + \ 99 - ((PHI_J_K[j][bval] + c) & mod_circ) 100 - chk = row_idx * M_SUB + c 101 - var = col_idx * M_SUB + pi 102 - edges[chk].append(var) 103 - 104 - process_layer(PROTO_0) 105 - process_layer(PROTO_1) 106 - process_layer(PROTO_2) 107 - 108 - # Sort and deduplicate 109 - h_rows = [sorted(set(row)) for row in edges] 110 - return h_rows 111 - 112 - 113 - def systematic_form(h_rows, m, n): 114 - """Gaussian elimination over GF(2) to get H into [A | I_m] form. 115 - 116 - Returns (perm, gen_p) where: 117 - - perm[i] = original column at position i 118 - - gen_p[j] = list of data column indices for parity bit j 119 - """ 120 - k = n - m 121 - 122 - # Build dense matrix 123 - mat = np.zeros((m, n), dtype=np.int8) 124 - for i, cols in enumerate(h_rows): 125 - for j in cols: 126 - mat[i, j] = 1 127 - 128 - perm = list(range(n)) 129 - 130 - def swap_cols(c1, c2): 131 - if c1 != c2: 132 - perm[c1], perm[c2] = perm[c2], perm[c1] 133 - mat[:, [c1, c2]] = mat[:, [c2, c1]] 134 - 135 - # Put identity in the last m columns 136 - for pivot in range(m): 137 - col = k + pivot 138 - # Find a row >= pivot with a 1 in column col 139 - found_row = -1 140 - for r in range(pivot, m): 141 - if mat[r, col] == 1: 142 - found_row = r 143 - break 144 - 145 - if found_row < 0: 146 - # Find a column in 0..k-1 with a 1 in some row >= pivot 147 - for c in range(k): 148 - for r in range(pivot, m): 149 - if mat[r, c] == 1: 150 - swap_cols(c, col) 151 - found_row = r 152 - break 153 - if found_row >= 0: 154 - break 155 - 156 - if found_row >= 0: 157 - if found_row != pivot: 158 - mat[[pivot, found_row]] = mat[[found_row, pivot]] 159 - 160 - # Eliminate 161 - for r2 in range(m): 162 - if r2 != pivot and mat[r2, col] == 1: 163 - mat[r2] ^= mat[pivot] 164 - 165 - # Extract generator parity sub-matrix 166 - gen_p = [] 167 - for j in range(m): 168 - cols = [i for i in range(k) if mat[j, i] == 1] 169 - gen_p.append(cols) 170 - 171 - return perm, gen_p 172 - 173 - 174 - def ldpc_encode(data_bytes, perm, gen_p): 175 - """Systematic LDPC encode. 176 - 177 - Args: 178 - data_bytes: 128 bytes (k=1024 bits) 179 - perm: column permutation from systematic_form 180 - gen_p: generator parity matrix from systematic_form 181 - 182 - Returns: 183 - 256 bytes (n=2048 transmitted bits) 184 - """ 185 - assert len(data_bytes) == K // 8, \ 186 - f"Expected {K // 8} bytes, got {len(data_bytes)}" 187 - 188 - # Extract k information bits (MSB-first) 189 - info = np.zeros(K, dtype=int) 190 - for i in range(K): 191 - byte_idx = i // 8 192 - bit_pos = 7 - (i % 8) 193 - info[i] = (data_bytes[byte_idx] >> bit_pos) & 1 194 - 195 - # Compute parity bits 196 - parity = np.zeros(H_M, dtype=int) 197 - for j in range(H_M): 198 - s = 0 199 - for i in gen_p[j]: 200 - s ^= info[i] 201 - parity[j] = s 202 - 203 - # Build full codeword [info | parity] 204 - full_bits = np.concatenate([info, parity]) 205 - 206 - # Puncture: keep only first N bits 207 - tx_bits = full_bits[:N] 208 - 209 - # Pack into bytes (MSB-first) 210 - result = [] 211 - for i in range(0, N, 8): 212 - byte = 0 213 - for j in range(8): 214 - if i + j < N: 215 - byte |= int(tx_bits[i + j]) << (7 - j) 216 - result.append(byte) 217 - return bytes(result) 218 - 219 - 220 - # --- Generate test vectors --- 221 - 222 - TRACE_DIR = (sys.argv[1] if len(sys.argv) > 1 223 - else os.path.join(os.path.dirname(__file__), "..", "traces")) 224 - 225 - # Test input patterns 226 - DATA_BYTES_LEN = K // 8 # 128 bytes 227 - 228 - inputs = [ 229 - ("zeros", bytes(DATA_BYTES_LEN)), 230 - ("ones", bytes([0xFF] * DATA_BYTES_LEN)), 231 - ("alternating_01", bytes([0x55] * DATA_BYTES_LEN)), 232 - ("alternating_10", bytes([0xAA] * DATA_BYTES_LEN)), 233 - ("ascending", bytes([i & 0xFF for i in range(DATA_BYTES_LEN)])), 234 - ] 235 - 236 - if __name__ == '__main__': 237 - print("Building CCSDS AR4JA H-matrix...", file=sys.stderr) 238 - h_rows = build_ar4ja_h() 239 - print("Computing systematic form...", file=sys.stderr) 240 - perm, gen_p = systematic_form(h_rows, H_M, H_N) 241 - print("Generating test vectors...", file=sys.stderr) 242 - 243 - os.makedirs(TRACE_DIR, exist_ok=True) 244 - with open(os.path.join(TRACE_DIR, "vectors.csv"), "w") as f: 245 - f.write("name,input_hex,encoded_hex\n") 246 - for name, data in inputs: 247 - coded = ldpc_encode(data, perm, gen_p) 248 - f.write(f"{name},{data.hex()},{coded.hex()}\n") 249 - 250 - print(f"Wrote {TRACE_DIR}/vectors.csv ({len(inputs)} vectors)") 251 - 252 - # Verify systematic property: first 128 bytes of codeword = data 253 - for name, data in inputs: 254 - coded = ldpc_encode(data, perm, gen_p) 255 - if coded[:128] != data: 256 - print(f"WARNING: {name}: systematic property violated!", file=sys.stderr) 257 - else: 258 - print(f" {name}: systematic property OK", file=sys.stderr)
-13
test/interop/commpy/scripts/generate.sh
··· 1 - #!/bin/bash 2 - set -euo pipefail 3 - SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)" 4 - TRACE_DIR="$SCRIPT_DIR/../traces" 5 - mkdir -p "$TRACE_DIR" 6 - TRACE_DIR="$(cd "$TRACE_DIR" && pwd)" 7 - 8 - cd "$SCRIPT_DIR" 9 - if [ ! -d .venv ]; then 10 - python3 -m venv .venv 11 - .venv/bin/pip install -r requirements.txt 12 - fi 13 - .venv/bin/python3 generate.py "$TRACE_DIR"
-3
test/interop/commpy/scripts/requirements.txt
··· 1 - scikit-commpy==0.8.0 2 - numpy 3 - scipy
-87
test/interop/commpy/test.ml
··· 1 - (** CCSDS LDPC interop tests for ocaml-ldpc. 2 - 3 - Traces generated by: Python CCSDS AR4JA encoder (same H-matrix 4 - construction as the OCaml implementation) 5 - Regenerate: dune build @regen-traces *) 6 - 7 - let trace path = Filename.concat "traces" path 8 - 9 - let bytes_of_hex hex = 10 - let len = String.length hex / 2 in 11 - Bytes.init len (fun i -> 12 - Char.chr (int_of_string ("0x" ^ String.sub hex (i * 2) 2))) 13 - 14 - let hex_of_bytes b = 15 - let buf = Buffer.create (Bytes.length b * 2) in 16 - Bytes.iter 17 - (fun c -> Buffer.add_string buf (Printf.sprintf "%02x" (Char.code c))) 18 - b; 19 - Buffer.contents buf 20 - 21 - type vector = { name : string; input : bytes; expected_coded : bytes } 22 - 23 - let vector_codec = 24 - Csvt.( 25 - Row.( 26 - obj (fun name input_hex encoded_hex -> 27 - { 28 - name; 29 - input = bytes_of_hex input_hex; 30 - expected_coded = bytes_of_hex encoded_hex; 31 - }) 32 - |> col "name" string ~enc:(fun v -> v.name) 33 - |> col "input_hex" string ~enc:(fun v -> hex_of_bytes v.input) 34 - |> col "encoded_hex" string ~enc:(fun v -> hex_of_bytes v.expected_coded) 35 - |> finish)) 36 - 37 - let parse_vectors () = 38 - match Csvt.decode_file vector_codec (trace "vectors.csv") with 39 - | Ok rows -> rows 40 - | Error e -> Alcotest.failf "CSV: %a" Csvt.pp_error e 41 - 42 - let ldpc = Ldpc.ccsds_rate_1_2 43 - 44 - let test_encode vec () = 45 - let coded = Ldpc.encode ldpc vec.input in 46 - if coded <> vec.expected_coded then 47 - Alcotest.failf "%s: encode mismatch\n expected: %s\n got: %s" 48 - vec.name 49 - (hex_of_bytes vec.expected_coded) 50 - (hex_of_bytes coded) 51 - 52 - let test_roundtrip vec () = 53 - let coded = Ldpc.encode ldpc vec.input in 54 - match Ldpc.decode ldpc coded with 55 - | Ok recovered -> 56 - if recovered <> vec.input then 57 - Alcotest.failf "%s: roundtrip mismatch\n expected: %s\n got: %s" 58 - vec.name (hex_of_bytes vec.input) (hex_of_bytes recovered) 59 - | Error e -> Alcotest.failf "%s: decode failed: %s" vec.name e 60 - 61 - let test_systematic vec () = 62 - let coded = Ldpc.encode ldpc vec.input in 63 - (* For the rate 1/2 code: first k/8 = 128 bytes should be the data *) 64 - let k_bytes = 128 in 65 - let prefix = Bytes.sub coded 0 k_bytes in 66 - if prefix <> vec.input then 67 - Alcotest.failf 68 - "%s: systematic property violated\n data: %s\n prefix: %s" vec.name 69 - (hex_of_bytes vec.input) (hex_of_bytes prefix) 70 - 71 - let () = 72 - let vectors = parse_vectors () in 73 - Alcotest.run "ldpc-interop" 74 - [ 75 - ( "encode", 76 - List.map 77 - (fun v -> Alcotest.test_case v.name `Quick (test_encode v)) 78 - vectors ); 79 - ( "systematic", 80 - List.map 81 - (fun v -> Alcotest.test_case v.name `Quick (test_systematic v)) 82 - vectors ); 83 - ( "roundtrip", 84 - List.map 85 - (fun v -> Alcotest.test_case v.name `Quick (test_roundtrip v)) 86 - vectors ); 87 - ]