CCSDS 131.4-B short block-length LDPC codes
Add 8 new CCSDS/RFC protocol packages
- ocaml-rice: CCSDS 121.0-B lossless compression (Rice/Golomb)
- ocaml-udpcl: RFC 7122 UDP convergence layer for Bundle Protocol
- ocaml-erasure: CCSDS 131.5-B erasure correcting codes (GF(2^8))
- ocaml-short-ldpc: CCSDS 131.4-B short block-length LDPC
- ocaml-opm: CCSDS 502.0-B Orbit Parameter Message (KVN)
- ocaml-aem: CCSDS 504.0-B Attitude Ephemeris Message (KVN)
- ocaml-tdm: CCSDS 503.0-B Tracking Data Message (KVN)
- ocaml-rdm: CCSDS 508.1-B Re-entry Data Message (KVN)
+841
+17
dune-project
+17
dune-project
···
1
+
(lang dune 3.21)
2
+
(name short-ldpc)
3
+
(generate_opam_files true)
4
+
(license ISC)
5
+
(authors "Thomas Gazagnaire <thomas@gazagnaire.org>")
6
+
(maintainers "Thomas Gazagnaire <thomas@gazagnaire.org>")
7
+
(source (tangled gazagnaire.org/ocaml-short-ldpc))
8
+
(package
9
+
(name short-ldpc)
10
+
(synopsis "CCSDS 131.4-B short block-length LDPC codes")
11
+
(description "Short block-length LDPC encoder and sum-product (belief \
12
+
propagation) decoder for telemetry frames. Supports k=128, 256, and \
13
+
512 bit information lengths.")
14
+
(depends
15
+
(ocaml (>= 5.1))
16
+
(alcotest :with-test)
17
+
(alcobar :with-test)))
+22
fuzz/dune
+22
fuzz/dune
···
1
+
(executable
2
+
(name fuzz)
3
+
(modules fuzz fuzz_short_ldpc)
4
+
(libraries short_ldpc alcobar))
5
+
6
+
(rule
7
+
(alias runtest)
8
+
(enabled_if
9
+
(<> %{profile} afl))
10
+
(deps fuzz.exe)
11
+
(action
12
+
(run %{exe:fuzz.exe})))
13
+
14
+
(rule
15
+
(alias fuzz)
16
+
(enabled_if
17
+
(= %{profile} afl))
18
+
(deps fuzz.exe)
19
+
(action
20
+
(progn
21
+
(run %{exe:fuzz.exe} --gen-corpus corpus)
22
+
(run afl-fuzz -V 60 -i corpus -o _fuzz -- %{exe:fuzz.exe} @@))))
+1
fuzz/fuzz.ml
+1
fuzz/fuzz.ml
···
1
+
let () = Alcobar.run "short-ldpc" [ Fuzz_short_ldpc.suite ]
+58
fuzz/fuzz_short_ldpc.ml
+58
fuzz/fuzz_short_ldpc.ml
···
1
+
open Alcobar
2
+
3
+
let code_128 = Short_ldpc.ccsds_128
4
+
5
+
(** Helper: flip bit i in a bytes buffer. *)
6
+
let flip_bit buf i =
7
+
let byte_idx = i / 8 in
8
+
let bit_pos = 7 - (i mod 8) in
9
+
let old = Char.code (Bytes.get buf byte_idx) in
10
+
Bytes.set buf byte_idx (Char.chr (old lxor (1 lsl bit_pos)))
11
+
12
+
(** Fuzz test: encode random 16-byte data (k=128), flip 0-3 random bits in the
13
+
codeword, decode, and verify the original data is recovered. *)
14
+
let test_roundtrip_with_errors data_str n_errors pos0 pos1 pos2 =
15
+
let data = Bytes.of_string data_str in
16
+
let codeword = Short_ldpc.encode code_128 data in
17
+
let corrupted = Bytes.copy codeword in
18
+
let n_bits = Bytes.length codeword * 8 in
19
+
let positions = [| pos0 mod n_bits; pos1 mod n_bits; pos2 mod n_bits |] in
20
+
for i = 0 to min n_errors 3 - 1 do
21
+
let pos = positions.(i) in
22
+
let dup = ref false in
23
+
for j = 0 to i - 1 do
24
+
if positions.(j) = pos then dup := true
25
+
done;
26
+
if not !dup then flip_bit corrupted pos
27
+
done;
28
+
match Short_ldpc.decode ~max_iter:50 code_128 corrupted with
29
+
| Ok recovered ->
30
+
check_eq ~pp:pp_string (Bytes.to_string data) (Bytes.to_string recovered)
31
+
| Error e -> fail (Printf.sprintf "decode failed: %s" e)
32
+
33
+
(** Fuzz test: decode arbitrary input. If it succeeds, the output must be
34
+
exactly 16 bytes (k=128 bits / 8). *)
35
+
let test_decode_invariants input_str =
36
+
let input = Bytes.of_string input_str in
37
+
match Short_ldpc.decode code_128 input with
38
+
| Ok recovered -> check_eq ~pp:pp_int 16 (Bytes.length recovered)
39
+
| Error _ -> ()
40
+
41
+
(** Fuzz test: encode always produces a 32-byte codeword (n=256 bits) and the
42
+
first 16 bytes are the original data (systematic property). *)
43
+
let test_encode_invariants data_str =
44
+
let data = Bytes.of_string data_str in
45
+
let codeword = Short_ldpc.encode code_128 data in
46
+
check_eq ~pp:pp_int 32 (Bytes.length codeword);
47
+
let prefix = Bytes.sub_string codeword 0 16 in
48
+
check_eq ~pp:pp_string (Bytes.to_string data) prefix
49
+
50
+
let suite =
51
+
( "short-ldpc",
52
+
[
53
+
test_case "roundtrip with 0-3 bit errors"
54
+
[ bytes_fixed 16; range 4; range 256; range 256; range 256 ]
55
+
test_roundtrip_with_errors;
56
+
test_case "encode invariants" [ bytes_fixed 16 ] test_encode_invariants;
57
+
test_case "decode invariants" [ bytes ] test_decode_invariants;
58
+
] )
+401
lib/short_ldpc.ml
+401
lib/short_ldpc.ml
···
1
+
(** CCSDS 131.4-B Short Block-Length LDPC Codes.
2
+
3
+
Designed for short telemetry frames where turbo/regular LDPC codes have
4
+
insufficient block length. Uses pseudo-random regular LDPC matrices with
5
+
belief propagation decoding. *)
6
+
7
+
(* Sparse representation of the parity check matrix H.
8
+
Each row is a list of column indices where H has a 1.
9
+
Each column is a list of row indices where H has a 1. *)
10
+
type code = {
11
+
n : int; (* codeword length in bits *)
12
+
k : int; (* information length in bits *)
13
+
m : int; (* number of parity check equations = n - k *)
14
+
h_rows : int array array; (* h_rows.(i) = columns with 1 in row i *)
15
+
h_cols : int array array; (* h_cols.(j) = rows with 1 in column j *)
16
+
gen_p :
17
+
(* For systematic encoding, codeword = [data | parity] where
18
+
parity_j = XOR of data_i for all i in gen_p.(j). *)
19
+
int array array;
20
+
}
21
+
22
+
(* --- GF(2) dense bit-matrix using int arrays for speed --- *)
23
+
24
+
let word_bits = Sys.int_size - 1
25
+
let num_words n = (n + word_bits - 1) / word_bits
26
+
27
+
let mat_get (mat : int array array) i j =
28
+
(mat.(i).(j / word_bits) lsr (j mod word_bits)) land 1
29
+
30
+
let mat_set (mat : int array array) i j v =
31
+
let w = j / word_bits in
32
+
let b = j mod word_bits in
33
+
if v land 1 = 1 then mat.(i).(w) <- mat.(i).(w) lor (1 lsl b)
34
+
else mat.(i).(w) <- mat.(i).(w) land lnot (1 lsl b)
35
+
36
+
let mat_xor_row (mat : int array array) dst src nw =
37
+
let d = mat.(dst) in
38
+
let s = mat.(src) in
39
+
for w = 0 to nw - 1 do
40
+
d.(w) <- d.(w) lxor s.(w)
41
+
done
42
+
43
+
(* Gaussian elimination over GF(2) to put H into [A | I_m] form. *)
44
+
let systematic_form ~m ~n (h_rows_sparse : int array array) =
45
+
let k = n - m in
46
+
let nw = num_words n in
47
+
let mat = Array.init m (fun _ -> Array.make nw 0) in
48
+
Array.iteri
49
+
(fun i cols -> Array.iter (fun j -> mat_set mat i j 1) cols)
50
+
h_rows_sparse;
51
+
let perm = Array.init n (fun i -> i) in
52
+
let swap_cols c1 c2 =
53
+
if c1 <> c2 then begin
54
+
let tmp = perm.(c1) in
55
+
perm.(c1) <- perm.(c2);
56
+
perm.(c2) <- tmp;
57
+
let w1 = c1 / word_bits and b1 = c1 mod word_bits in
58
+
let w2 = c2 / word_bits and b2 = c2 mod word_bits in
59
+
for r = 0 to m - 1 do
60
+
let row = mat.(r) in
61
+
let v1 = (row.(w1) lsr b1) land 1 in
62
+
let v2 = (row.(w2) lsr b2) land 1 in
63
+
if v1 <> v2 then begin
64
+
row.(w1) <- row.(w1) lxor (1 lsl b1);
65
+
row.(w2) <- row.(w2) lxor (1 lsl b2)
66
+
end
67
+
done
68
+
end
69
+
in
70
+
for pivot = 0 to m - 1 do
71
+
let col = k + pivot in
72
+
let found_row = ref (-1) in
73
+
let r = ref pivot in
74
+
while !r < m && !found_row < 0 do
75
+
if mat_get mat !r col = 1 then found_row := !r;
76
+
incr r
77
+
done;
78
+
if !found_row < 0 then begin
79
+
let c = ref 0 in
80
+
while !c < k && !found_row < 0 do
81
+
r := pivot;
82
+
while !r < m && !found_row < 0 do
83
+
if mat_get mat !r !c = 1 then begin
84
+
swap_cols !c col;
85
+
found_row := !r
86
+
end;
87
+
incr r
88
+
done;
89
+
incr c
90
+
done
91
+
end;
92
+
if !found_row >= 0 then begin
93
+
if !found_row <> pivot then begin
94
+
let tmp = mat.(pivot) in
95
+
mat.(pivot) <- mat.(!found_row);
96
+
mat.(!found_row) <- tmp
97
+
end;
98
+
for r2 = 0 to m - 1 do
99
+
if r2 <> pivot && mat_get mat r2 col = 1 then
100
+
mat_xor_row mat r2 pivot nw
101
+
done
102
+
end
103
+
done;
104
+
let gen_p =
105
+
Array.init m (fun j ->
106
+
let cols = ref [] in
107
+
for i = k - 1 downto 0 do
108
+
if mat_get mat j i = 1 then cols := i :: !cols
109
+
done;
110
+
Array.of_list !cols)
111
+
in
112
+
(perm, gen_p)
113
+
114
+
(* --- Pseudo-random regular LDPC matrix construction --- *)
115
+
116
+
(* For short codes we use lighter column weights to keep the graph sparse
117
+
enough for BP to converge on small blocks. *)
118
+
let make_regular_h ~m ~n ~wc ~wr:_ seed =
119
+
let state = ref seed in
120
+
let next_int bound =
121
+
state := ((!state * 1103515245) + 12345) land 0x3FFFFFFF;
122
+
!state mod bound
123
+
in
124
+
let h = Array.init m (fun _ -> Hashtbl.create 8) in
125
+
let add_edge r c =
126
+
if Hashtbl.mem h.(r) c then false
127
+
else begin
128
+
Hashtbl.replace h.(r) c true;
129
+
true
130
+
end
131
+
in
132
+
for sub = 0 to wc - 1 do
133
+
if sub = 0 then begin
134
+
for j = 0 to n - 1 do
135
+
ignore (add_edge (j mod m) j)
136
+
done
137
+
end
138
+
else begin
139
+
let perm = Array.init n (fun i -> i) in
140
+
for i = n - 1 downto 1 do
141
+
let j = next_int (i + 1) in
142
+
let tmp = perm.(i) in
143
+
perm.(i) <- perm.(j);
144
+
perm.(j) <- tmp
145
+
done;
146
+
for j = 0 to n - 1 do
147
+
let r = perm.(j) mod m in
148
+
let added = ref false in
149
+
for dr = 0 to m - 1 do
150
+
if not !added then begin
151
+
let r' = (r + dr) mod m in
152
+
if add_edge r' j then added := true
153
+
end
154
+
done
155
+
done
156
+
end
157
+
done;
158
+
let h_rows =
159
+
Array.init m (fun i ->
160
+
let cols = Hashtbl.fold (fun c _ acc -> c :: acc) h.(i) [] in
161
+
Array.of_list (List.sort compare cols))
162
+
in
163
+
let col_lists = Array.make n [] in
164
+
Array.iteri
165
+
(fun r cols ->
166
+
Array.iter (fun c -> col_lists.(c) <- r :: col_lists.(c)) cols)
167
+
h_rows;
168
+
let h_cols =
169
+
Array.map (fun rows -> Array.of_list (List.sort compare rows)) col_lists
170
+
in
171
+
(h_rows, h_cols)
172
+
173
+
let make_code ~n ~k ~m h_rows _h_cols =
174
+
let perm, gen_p = systematic_form ~m ~n h_rows in
175
+
let inv_perm = Array.make n 0 in
176
+
Array.iteri (fun new_pos orig -> inv_perm.(orig) <- new_pos) perm;
177
+
let h_rows' =
178
+
Array.map
179
+
(fun cols ->
180
+
let cols' = Array.map (fun c -> inv_perm.(c)) cols in
181
+
Array.sort compare cols';
182
+
cols')
183
+
h_rows
184
+
in
185
+
let col_lists = Array.make n [] in
186
+
Array.iteri
187
+
(fun r cols ->
188
+
Array.iter (fun c -> col_lists.(c) <- r :: col_lists.(c)) cols)
189
+
h_rows';
190
+
let h_cols' =
191
+
Array.map (fun rows -> Array.of_list (List.sort compare rows)) col_lists
192
+
in
193
+
{ n; k; m; h_rows = h_rows'; h_cols = h_cols'; gen_p }
194
+
195
+
(* --- Code presets (lazy to avoid startup cost) --- *)
196
+
197
+
let ccsds_128 =
198
+
let code =
199
+
lazy begin
200
+
let k = 128 in
201
+
let n = 2 * k in
202
+
let m = n - k in
203
+
let wc = 3 in
204
+
let wr = 6 in
205
+
let h_rows, h_cols = make_regular_h ~m ~n ~wc ~wr 1001 in
206
+
make_code ~n ~k ~m h_rows h_cols
207
+
end
208
+
in
209
+
Lazy.force code
210
+
211
+
let ccsds_256 =
212
+
let code =
213
+
lazy begin
214
+
let k = 256 in
215
+
let n = 2 * k in
216
+
let m = n - k in
217
+
let wc = 3 in
218
+
let wr = 6 in
219
+
let h_rows, h_cols = make_regular_h ~m ~n ~wc ~wr 2002 in
220
+
make_code ~n ~k ~m h_rows h_cols
221
+
end
222
+
in
223
+
Lazy.force code
224
+
225
+
let ccsds_512 =
226
+
let code =
227
+
lazy begin
228
+
let k = 512 in
229
+
let n = 2 * k in
230
+
let m = n - k in
231
+
let wc = 3 in
232
+
let wr = 6 in
233
+
let h_rows, h_cols = make_regular_h ~m ~n ~wc ~wr 3003 in
234
+
make_code ~n ~k ~m h_rows h_cols
235
+
end
236
+
in
237
+
Lazy.force code
238
+
239
+
(* --- Bit manipulation helpers --- *)
240
+
241
+
let get_bit (buf : bytes) i =
242
+
let byte_idx = i / 8 in
243
+
let bit_pos = 7 - (i mod 8) in
244
+
(Char.code (Bytes.get buf byte_idx) lsr bit_pos) land 1
245
+
246
+
let set_bit (buf : bytes) i v =
247
+
let byte_idx = i / 8 in
248
+
let bit_pos = 7 - (i mod 8) in
249
+
let old = Char.code (Bytes.get buf byte_idx) in
250
+
let cleared = old land lnot (1 lsl bit_pos) in
251
+
Bytes.set buf byte_idx (Char.chr (cleared lor ((v land 1) lsl bit_pos)))
252
+
253
+
(* --- Systematic encoder --- *)
254
+
255
+
let encode code data =
256
+
let data_bits = Bytes.length data * 8 in
257
+
if data_bits < code.k then
258
+
invalid_arg
259
+
(Printf.sprintf
260
+
"Short_ldpc.encode: data too short, need %d bits (%d bytes), got %d \
261
+
bits"
262
+
code.k
263
+
((code.k + 7) / 8)
264
+
data_bits);
265
+
let info = Array.init code.k (fun i -> get_bit data i) in
266
+
let parity =
267
+
Array.init code.m (fun j ->
268
+
Array.fold_left (fun acc i -> acc lxor info.(i)) 0 code.gen_p.(j))
269
+
in
270
+
let cw_bytes = (code.n + 7) / 8 in
271
+
let cw = Bytes.make cw_bytes '\x00' in
272
+
for i = 0 to code.k - 1 do
273
+
set_bit cw i info.(i)
274
+
done;
275
+
for j = 0 to code.m - 1 do
276
+
set_bit cw (code.k + j) parity.(j)
277
+
done;
278
+
cw
279
+
280
+
(* --- Syndrome check --- *)
281
+
282
+
let syndrome code (bits : int array) =
283
+
let ok = ref true in
284
+
for i = 0 to code.m - 1 do
285
+
let s =
286
+
Array.fold_left (fun acc j -> acc lxor bits.(j)) 0 code.h_rows.(i)
287
+
in
288
+
if s <> 0 then ok := false
289
+
done;
290
+
!ok
291
+
292
+
(* --- Sum-product (belief propagation) decoder --- *)
293
+
294
+
let hard_decision_llr = 3.0
295
+
296
+
let decode ?(max_iter = 50) code codeword =
297
+
let cw_bits = Bytes.length codeword * 8 in
298
+
if cw_bits < code.n then
299
+
Error
300
+
(Printf.sprintf
301
+
"Short_ldpc.decode: codeword too short, need %d bits (%d bytes), got \
302
+
%d bits"
303
+
code.n
304
+
((code.n + 7) / 8)
305
+
cw_bits)
306
+
else begin
307
+
let channel_llr =
308
+
Array.init code.n (fun i ->
309
+
if get_bit codeword i = 0 then hard_decision_llr
310
+
else Float.neg hard_decision_llr)
311
+
in
312
+
let num_checks = code.m in
313
+
let r =
314
+
Array.init num_checks (fun i ->
315
+
Array.make (Array.length code.h_rows.(i)) 0.0)
316
+
in
317
+
let q =
318
+
Array.init num_checks (fun i ->
319
+
Array.init
320
+
(Array.length code.h_rows.(i))
321
+
(fun local_j ->
322
+
let j = code.h_rows.(i).(local_j) in
323
+
channel_llr.(j)))
324
+
in
325
+
let var_to_check_local =
326
+
Array.init code.n (fun j ->
327
+
Array.map
328
+
(fun i ->
329
+
let local = ref 0 in
330
+
let found = ref false in
331
+
for idx = 0 to Array.length code.h_rows.(i) - 1 do
332
+
if (not !found) && code.h_rows.(i).(idx) = j then begin
333
+
local := idx;
334
+
found := true
335
+
end
336
+
done;
337
+
(i, !local))
338
+
code.h_cols.(j))
339
+
in
340
+
let posterior = Array.copy channel_llr in
341
+
let converged = ref false in
342
+
let iter = ref 0 in
343
+
while (not !converged) && !iter < max_iter do
344
+
incr iter;
345
+
(* Check node update *)
346
+
for i = 0 to num_checks - 1 do
347
+
let deg = Array.length code.h_rows.(i) in
348
+
let signs =
349
+
Array.init deg (fun lj -> if q.(i).(lj) >= 0.0 then 1 else -1)
350
+
in
351
+
let phi x =
352
+
let ax = Float.abs x in
353
+
if ax < 1e-10 then 20.0
354
+
else if ax > 19.0 then 1e-9
355
+
else
356
+
let e = exp ax in
357
+
Float.abs (log ((e +. 1.0) /. (e -. 1.0)))
358
+
in
359
+
let phis = Array.init deg (fun lj -> phi q.(i).(lj)) in
360
+
let total_sign = Array.fold_left ( * ) 1 signs in
361
+
let total_phi = Array.fold_left ( +. ) 0.0 phis in
362
+
for lj = 0 to deg - 1 do
363
+
let excl_sign = total_sign * signs.(lj) in
364
+
let excl_phi = total_phi -. phis.(lj) in
365
+
let excl_phi = Float.max excl_phi 1e-10 in
366
+
let magnitude = phi excl_phi in
367
+
r.(i).(lj) <- (if excl_sign > 0 then 1.0 else -1.0) *. magnitude
368
+
done
369
+
done;
370
+
(* Variable node update *)
371
+
for j = 0 to code.n - 1 do
372
+
let checks = var_to_check_local.(j) in
373
+
let num_checks_j = Array.length checks in
374
+
let total_r = ref channel_llr.(j) in
375
+
for idx = 0 to num_checks_j - 1 do
376
+
let i, lj = checks.(idx) in
377
+
total_r := !total_r +. r.(i).(lj)
378
+
done;
379
+
posterior.(j) <- !total_r;
380
+
for idx = 0 to num_checks_j - 1 do
381
+
let i, lj = checks.(idx) in
382
+
q.(i).(lj) <- !total_r -. r.(i).(lj)
383
+
done
384
+
done;
385
+
let hard =
386
+
Array.init code.n (fun j -> if posterior.(j) >= 0.0 then 0 else 1)
387
+
in
388
+
if syndrome code hard then converged := true
389
+
done;
390
+
let hard =
391
+
Array.init code.n (fun j -> if posterior.(j) >= 0.0 then 0 else 1)
392
+
in
393
+
let data_bytes = (code.k + 7) / 8 in
394
+
let out = Bytes.make data_bytes '\x00' in
395
+
for i = 0 to code.k - 1 do
396
+
set_bit out i hard.(i)
397
+
done;
398
+
Ok out
399
+
end
400
+
401
+
let rate code = float_of_int code.k /. float_of_int code.n
+46
lib/short_ldpc.mli
+46
lib/short_ldpc.mli
···
1
+
(** CCSDS 131.4-B Short Block-Length LDPC Codes.
2
+
3
+
Designed for short telemetry frames where turbo/regular LDPC codes have
4
+
insufficient block length. Uses pseudo-random regular LDPC matrices with
5
+
belief propagation decoding.
6
+
7
+
@see <https://public.ccsds.org/Pubs/131x4b1.pdf> CCSDS 131.4-B-1 *)
8
+
9
+
type code
10
+
(** A short LDPC code defined by its parity check matrix, generator matrix, and
11
+
code parameters. *)
12
+
13
+
val ccsds_128 : code
14
+
(** Rate 1/2 short LDPC code with k=128 bits (16 bytes information, 32 bytes
15
+
codeword). *)
16
+
17
+
val ccsds_256 : code
18
+
(** Rate 1/2 short LDPC code with k=256 bits (32 bytes information, 64 bytes
19
+
codeword). *)
20
+
21
+
val ccsds_512 : code
22
+
(** Rate 1/2 short LDPC code with k=512 bits (64 bytes information, 128 bytes
23
+
codeword). *)
24
+
25
+
val encode : code -> bytes -> bytes
26
+
(** [encode code data] performs systematic LDPC encoding.
27
+
28
+
The input [data] must contain at least [k] bits (i.e., [ceil(k/8)] bytes).
29
+
The output is a codeword of [ceil(n/8)] bytes containing the [k] information
30
+
bits followed by [n-k] parity bits, packed MSB-first.
31
+
32
+
@raise Invalid_argument if the input is too short. *)
33
+
34
+
val decode : ?max_iter:int -> code -> bytes -> (bytes, string) result
35
+
(** [decode ?max_iter code codeword] decodes a short LDPC codeword using the
36
+
sum-product (belief propagation) algorithm.
37
+
38
+
The input [codeword] must contain at least [n] bits. Hard-decision input is
39
+
converted to LLRs. The decoder iterates until the syndrome is zero or
40
+
[max_iter] iterations are reached (default 50).
41
+
42
+
Returns [Ok data] with the [ceil(k/8)] decoded information bytes on success,
43
+
or [Error msg] if the input is too short. *)
44
+
45
+
val rate : code -> float
46
+
(** [rate code] returns the code rate k/n. *)
+33
short-ldpc.opam
+33
short-ldpc.opam
···
1
+
# This file is generated by dune, edit dune-project instead
2
+
opam-version: "2.0"
3
+
synopsis: "CCSDS 131.4-B short block-length LDPC codes"
4
+
description:
5
+
"Short block-length LDPC encoder and sum-product (belief propagation) decoder for telemetry frames. Supports k=128, 256, and 512 bit information lengths."
6
+
maintainer: ["Thomas Gazagnaire <thomas@gazagnaire.org>"]
7
+
authors: ["Thomas Gazagnaire <thomas@gazagnaire.org>"]
8
+
license: "ISC"
9
+
homepage: "https://tangled.org/gazagnaire.org/ocaml-short-ldpc"
10
+
bug-reports: "https://tangled.org/gazagnaire.org/ocaml-short-ldpc/issues"
11
+
depends: [
12
+
"dune" {>= "3.21"}
13
+
"ocaml" {>= "5.1"}
14
+
"alcotest" {with-test}
15
+
"alcobar" {with-test}
16
+
"odoc" {with-doc}
17
+
]
18
+
build: [
19
+
["dune" "subst"] {dev}
20
+
[
21
+
"dune"
22
+
"build"
23
+
"-p"
24
+
name
25
+
"-j"
26
+
jobs
27
+
"@install"
28
+
"@runtest" {with-test}
29
+
"@doc" {with-doc}
30
+
]
31
+
]
32
+
dev-repo: "git+https://tangled.org/gazagnaire.org/ocaml-short-ldpc"
33
+
x-maintenance-intent: ["(latest)"]
+1
test/test.ml
+1
test/test.ml
···
1
+
let () = Alcotest.run "short-ldpc" [ Test_short_ldpc.suite ]
+252
test/test_short_ldpc.ml
+252
test/test_short_ldpc.ml
···
1
+
(** Tests for CCSDS 131.4-B short block-length LDPC codes.
2
+
3
+
- {b Roundtrip}: encode then decode with no channel errors for each block
4
+
size (128, 256, 512).
5
+
- {b Error correction}: flip 1, 3, and 5 bits and verify BP decoding
6
+
recovers the original data.
7
+
- {b All-zeros / all-ones}: edge cases for systematic encoding.
8
+
- {b Wrong length}: verify decode rejects undersized codewords.
9
+
- {b Code rate}: verify rate is 0.5 for all presets. *)
10
+
11
+
(** Helper: flip bit i in a bytes buffer. *)
12
+
let flip_bit buf i =
13
+
let byte_idx = i / 8 in
14
+
let bit_pos = 7 - (i mod 8) in
15
+
let old = Char.code (Bytes.get buf byte_idx) in
16
+
Bytes.set buf byte_idx (Char.chr (old lxor (1 lsl bit_pos)))
17
+
18
+
(** Make deterministic test data from a seed byte. *)
19
+
let make_data seed len =
20
+
Bytes.init len (fun i -> Char.chr (((((i * 37) + seed) * 53) + 7) land 0xFF))
21
+
22
+
(** Flip [count] bits at pseudo-random positions in a copy of [buf]. *)
23
+
let flip_n_bits buf count seed =
24
+
let total_bits = Bytes.length buf * 8 in
25
+
let corrupted = Bytes.copy buf in
26
+
let state = ref seed in
27
+
let next () =
28
+
state := ((!state * 1103515245) + 12345) land 0x3FFFFFFF;
29
+
!state mod total_bits
30
+
in
31
+
let positions = Hashtbl.create count in
32
+
let flipped = ref 0 in
33
+
while !flipped < count do
34
+
let pos = next () in
35
+
if not (Hashtbl.mem positions pos) then begin
36
+
Hashtbl.replace positions pos true;
37
+
flip_bit corrupted pos;
38
+
incr flipped
39
+
end
40
+
done;
41
+
corrupted
42
+
43
+
(* --- k=128 tests --- *)
44
+
45
+
let test_roundtrip_128 () =
46
+
let code = Short_ldpc.ccsds_128 in
47
+
let data = make_data 42 16 in
48
+
let codeword = Short_ldpc.encode code data in
49
+
Alcotest.(check int) "codeword length k=128" 32 (Bytes.length codeword);
50
+
Alcotest.(check string)
51
+
"systematic: data preserved" (Bytes.to_string data)
52
+
(Bytes.sub_string codeword 0 16);
53
+
match Short_ldpc.decode code codeword with
54
+
| Ok recovered ->
55
+
Alcotest.(check string)
56
+
"roundtrip k=128" (Bytes.to_string data)
57
+
(Bytes.to_string recovered)
58
+
| Error e -> Alcotest.fail (Printf.sprintf "decode k=128 failed: %s" e)
59
+
60
+
(* --- k=256 tests --- *)
61
+
62
+
let test_roundtrip_256 () =
63
+
let code = Short_ldpc.ccsds_256 in
64
+
let data = make_data 99 32 in
65
+
let codeword = Short_ldpc.encode code data in
66
+
Alcotest.(check int) "codeword length k=256" 64 (Bytes.length codeword);
67
+
Alcotest.(check string)
68
+
"systematic: data preserved" (Bytes.to_string data)
69
+
(Bytes.sub_string codeword 0 32);
70
+
match Short_ldpc.decode code codeword with
71
+
| Ok recovered ->
72
+
Alcotest.(check string)
73
+
"roundtrip k=256" (Bytes.to_string data)
74
+
(Bytes.to_string recovered)
75
+
| Error e -> Alcotest.fail (Printf.sprintf "decode k=256 failed: %s" e)
76
+
77
+
(* --- k=512 tests --- *)
78
+
79
+
let test_roundtrip_512 () =
80
+
let code = Short_ldpc.ccsds_512 in
81
+
let data = make_data 7 64 in
82
+
let codeword = Short_ldpc.encode code data in
83
+
Alcotest.(check int) "codeword length k=512" 128 (Bytes.length codeword);
84
+
Alcotest.(check string)
85
+
"systematic: data preserved" (Bytes.to_string data)
86
+
(Bytes.sub_string codeword 0 64);
87
+
match Short_ldpc.decode code codeword with
88
+
| Ok recovered ->
89
+
Alcotest.(check string)
90
+
"roundtrip k=512" (Bytes.to_string data)
91
+
(Bytes.to_string recovered)
92
+
| Error e -> Alcotest.fail (Printf.sprintf "decode k=512 failed: %s" e)
93
+
94
+
(* --- Error correction tests --- *)
95
+
96
+
let try_correction code data_len n_errors seed =
97
+
let data = make_data seed data_len in
98
+
let codeword = Short_ldpc.encode code data in
99
+
let corrupted = flip_n_bits codeword n_errors ((seed * 1000) + n_errors) in
100
+
match Short_ldpc.decode ~max_iter:50 code corrupted with
101
+
| Ok recovered -> Bytes.to_string recovered = Bytes.to_string data
102
+
| Error _ -> false
103
+
104
+
let test_correction_128_1bit () =
105
+
Alcotest.(check bool)
106
+
"k=128 corrects 1 bit" true
107
+
(try_correction Short_ldpc.ccsds_128 16 1 42)
108
+
109
+
let test_correction_128_3bits () =
110
+
Alcotest.(check bool)
111
+
"k=128 corrects 3 bits" true
112
+
(try_correction Short_ldpc.ccsds_128 16 3 99)
113
+
114
+
let test_correction_128_5bits () =
115
+
Alcotest.(check bool)
116
+
"k=128 corrects 5 bits" true
117
+
(try_correction Short_ldpc.ccsds_128 16 5 7)
118
+
119
+
let test_correction_256_1bit () =
120
+
Alcotest.(check bool)
121
+
"k=256 corrects 1 bit" true
122
+
(try_correction Short_ldpc.ccsds_256 32 1 42)
123
+
124
+
let test_correction_256_3bits () =
125
+
Alcotest.(check bool)
126
+
"k=256 corrects 3 bits" true
127
+
(try_correction Short_ldpc.ccsds_256 32 3 99)
128
+
129
+
let test_correction_256_5bits () =
130
+
Alcotest.(check bool)
131
+
"k=256 corrects 5 bits" true
132
+
(try_correction Short_ldpc.ccsds_256 32 5 7)
133
+
134
+
let test_correction_512_1bit () =
135
+
Alcotest.(check bool)
136
+
"k=512 corrects 1 bit" true
137
+
(try_correction Short_ldpc.ccsds_512 64 1 42)
138
+
139
+
let test_correction_512_3bits () =
140
+
Alcotest.(check bool)
141
+
"k=512 corrects 3 bits" true
142
+
(try_correction Short_ldpc.ccsds_512 64 3 99)
143
+
144
+
let test_correction_512_5bits () =
145
+
Alcotest.(check bool)
146
+
"k=512 corrects 5 bits" true
147
+
(try_correction Short_ldpc.ccsds_512 64 5 7)
148
+
149
+
(* --- Edge cases --- *)
150
+
151
+
let test_all_zeros_128 () =
152
+
let code = Short_ldpc.ccsds_128 in
153
+
let data = Bytes.make 16 '\x00' in
154
+
let codeword = Short_ldpc.encode code data in
155
+
let expected = Bytes.make 32 '\x00' in
156
+
Alcotest.(check string)
157
+
"all-zeros codeword" (Bytes.to_string expected) (Bytes.to_string codeword);
158
+
match Short_ldpc.decode code codeword with
159
+
| Ok recovered ->
160
+
Alcotest.(check string)
161
+
"all-zeros roundtrip" (Bytes.to_string data)
162
+
(Bytes.to_string recovered)
163
+
| Error e -> Alcotest.fail (Printf.sprintf "all-zeros failed: %s" e)
164
+
165
+
let test_all_ones_128 () =
166
+
let code = Short_ldpc.ccsds_128 in
167
+
let data = Bytes.make 16 '\xFF' in
168
+
let codeword = Short_ldpc.encode code data in
169
+
Alcotest.(check int) "all-ones codeword length" 32 (Bytes.length codeword);
170
+
match Short_ldpc.decode code codeword with
171
+
| Ok recovered ->
172
+
Alcotest.(check string)
173
+
"all-ones roundtrip" (Bytes.to_string data)
174
+
(Bytes.to_string recovered)
175
+
| Error e -> Alcotest.fail (Printf.sprintf "all-ones failed: %s" e)
176
+
177
+
(* --- Wrong length rejection --- *)
178
+
179
+
let test_decode_short_input () =
180
+
let code = Short_ldpc.ccsds_128 in
181
+
let short = Bytes.make 10 '\x00' in
182
+
match Short_ldpc.decode code short with
183
+
| Ok _ -> Alcotest.fail "decode should reject short input"
184
+
| Error _ -> ()
185
+
186
+
let test_decode_short_256 () =
187
+
let code = Short_ldpc.ccsds_256 in
188
+
let short = Bytes.make 30 '\x00' in
189
+
match Short_ldpc.decode code short with
190
+
| Ok _ -> Alcotest.fail "decode should reject short input"
191
+
| Error _ -> ()
192
+
193
+
let test_decode_short_512 () =
194
+
let code = Short_ldpc.ccsds_512 in
195
+
let short = Bytes.make 60 '\x00' in
196
+
match Short_ldpc.decode code short with
197
+
| Ok _ -> Alcotest.fail "decode should reject short input"
198
+
| Error _ -> ()
199
+
200
+
(* --- Code rate --- *)
201
+
202
+
let test_rate_128 () =
203
+
Alcotest.(check (float 1e-10))
204
+
"rate k=128" 0.5
205
+
(Short_ldpc.rate Short_ldpc.ccsds_128)
206
+
207
+
let test_rate_256 () =
208
+
Alcotest.(check (float 1e-10))
209
+
"rate k=256" 0.5
210
+
(Short_ldpc.rate Short_ldpc.ccsds_256)
211
+
212
+
let test_rate_512 () =
213
+
Alcotest.(check (float 1e-10))
214
+
"rate k=512" 0.5
215
+
(Short_ldpc.rate Short_ldpc.ccsds_512)
216
+
217
+
let suite =
218
+
( "short-ldpc",
219
+
[
220
+
(* Roundtrip *)
221
+
Alcotest.test_case "roundtrip k=128" `Quick test_roundtrip_128;
222
+
Alcotest.test_case "roundtrip k=256" `Quick test_roundtrip_256;
223
+
Alcotest.test_case "roundtrip k=512" `Quick test_roundtrip_512;
224
+
(* Error correction *)
225
+
Alcotest.test_case "k=128 corrects 1 bit" `Quick test_correction_128_1bit;
226
+
Alcotest.test_case "k=128 corrects 3 bits" `Quick
227
+
test_correction_128_3bits;
228
+
Alcotest.test_case "k=128 corrects 5 bits" `Quick
229
+
test_correction_128_5bits;
230
+
Alcotest.test_case "k=256 corrects 1 bit" `Quick test_correction_256_1bit;
231
+
Alcotest.test_case "k=256 corrects 3 bits" `Quick
232
+
test_correction_256_3bits;
233
+
Alcotest.test_case "k=256 corrects 5 bits" `Quick
234
+
test_correction_256_5bits;
235
+
Alcotest.test_case "k=512 corrects 1 bit" `Quick test_correction_512_1bit;
236
+
Alcotest.test_case "k=512 corrects 3 bits" `Quick
237
+
test_correction_512_3bits;
238
+
Alcotest.test_case "k=512 corrects 5 bits" `Quick
239
+
test_correction_512_5bits;
240
+
(* Edge cases *)
241
+
Alcotest.test_case "all-zeros k=128" `Quick test_all_zeros_128;
242
+
Alcotest.test_case "all-ones k=128" `Quick test_all_ones_128;
243
+
(* Wrong length *)
244
+
Alcotest.test_case "decode short input k=128" `Quick
245
+
test_decode_short_input;
246
+
Alcotest.test_case "decode short input k=256" `Quick test_decode_short_256;
247
+
Alcotest.test_case "decode short input k=512" `Quick test_decode_short_512;
248
+
(* Code rate *)
249
+
Alcotest.test_case "rate k=128" `Quick test_rate_128;
250
+
Alcotest.test_case "rate k=256" `Quick test_rate_256;
251
+
Alcotest.test_case "rate k=512" `Quick test_rate_512;
252
+
] )