LDPC codes with belief propagation decoding
Split generic coding algorithms into reusable packages
Extract from ocaml-tm-sync into standalone packages:
- ocaml-viterbi: convolutional codes + Viterbi decoder (6 tests)
- ocaml-turbo: parallel concatenated codes + BCJR/MAP decoder (10 tests)
- ocaml-ldpc: belief propagation decoder (7 tests)
These implement standard algorithms reusable beyond CCSDS:
Viterbi (DVB-T, DAB, 802.11a/g), turbo (3G, LTE),
LDPC (DVB-S2, 5G NR, WiFi 6).
tm-sync keeps CCSDS-specific code (randomizer, RS interleaving)
and re-exports the generic packages for backward compat.
+640
+1
.ocamlformat
+1
.ocamlformat
···
1
+
version = 0.28.1
+16
dune-project
+16
dune-project
···
1
+
(lang dune 3.21)
2
+
(name 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-ldpc))
8
+
(package
9
+
(name ldpc)
10
+
(synopsis "LDPC codes with belief propagation decoding")
11
+
(description "Low-Density Parity-Check encoder and sum-product (belief \
12
+
propagation) decoder. Configurable parity check matrix. \
13
+
Includes a CCSDS-style rate 1/2 preset for testing.")
14
+
(depends
15
+
(ocaml (>= 5.1))
16
+
(alcotest :with-test)))
+32
ldpc.opam
+32
ldpc.opam
···
1
+
# This file is generated by dune, edit dune-project instead
2
+
opam-version: "2.0"
3
+
synopsis: "LDPC codes with belief propagation decoding"
4
+
description:
5
+
"Low-Density Parity-Check encoder and sum-product (belief propagation) decoder. Configurable parity check matrix. Includes a CCSDS-style rate 1/2 preset for testing."
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-ldpc"
10
+
bug-reports: "https://tangled.org/gazagnaire.org/ocaml-ldpc/issues"
11
+
depends: [
12
+
"dune" {>= "3.21"}
13
+
"ocaml" {>= "5.1"}
14
+
"alcotest" {with-test}
15
+
"odoc" {with-doc}
16
+
]
17
+
build: [
18
+
["dune" "subst"] {dev}
19
+
[
20
+
"dune"
21
+
"build"
22
+
"-p"
23
+
name
24
+
"-j"
25
+
jobs
26
+
"@install"
27
+
"@runtest" {with-test}
28
+
"@doc" {with-doc}
29
+
]
30
+
]
31
+
dev-repo: "git+https://tangled.org/gazagnaire.org/ocaml-ldpc"
32
+
x-maintenance-intent: ["(latest)"]
+423
lib/ldpc.ml
+423
lib/ldpc.ml
···
1
+
(** Low-Density Parity-Check codes with belief propagation decoding.
2
+
3
+
Generic LDPC encoder/decoder framework with systematic encoding and
4
+
sum-product (belief propagation) decoding. *)
5
+
6
+
(* Sparse representation of the parity check matrix H.
7
+
Each row is a list of column indices where H has a 1.
8
+
Each column is a list of row indices where H has a 1. *)
9
+
type code = {
10
+
n : int; (* codeword length *)
11
+
k : int; (* information length *)
12
+
m : int; (* number of parity check equations = n - k *)
13
+
h_rows : int array array; (* h_rows.(i) = columns with 1 in row i *)
14
+
h_cols : int array array; (* h_cols.(j) = rows with 1 in column j *)
15
+
gen_p :
16
+
(* For systematic encoding, codeword = [data | parity] where
17
+
parity_j = XOR of data_i for all i in gen_p.(j). *)
18
+
int array array;
19
+
}
20
+
21
+
(* --- GF(2) dense bit-matrix using int arrays for speed --- *)
22
+
23
+
(* Each row is an array of ints, each int holds Sys.int_size - 1 bits.
24
+
We use 62-bit words on 64-bit OCaml. *)
25
+
let word_bits = Sys.int_size - 1 (* 62 or 30 *)
26
+
let num_words n = (n + word_bits - 1) / word_bits
27
+
28
+
let mat_get (mat : int array array) i j =
29
+
(mat.(i).(j / word_bits) lsr (j mod word_bits)) land 1
30
+
31
+
let mat_set (mat : int array array) i j v =
32
+
let w = j / word_bits in
33
+
let b = j mod word_bits in
34
+
if v land 1 = 1 then mat.(i).(w) <- mat.(i).(w) lor (1 lsl b)
35
+
else mat.(i).(w) <- mat.(i).(w) land lnot (1 lsl b)
36
+
37
+
let mat_xor_row (mat : int array array) dst src nw =
38
+
let d = mat.(dst) in
39
+
let s = mat.(src) in
40
+
for w = 0 to nw - 1 do
41
+
d.(w) <- d.(w) lxor s.(w)
42
+
done
43
+
44
+
(* Gaussian elimination over GF(2) to put H into [A | I_m] form.
45
+
We track column permutations without physically swapping bits --
46
+
instead we swap column indices in a virtual permutation array and
47
+
physically swap columns only in word-level bulk. *)
48
+
let systematic_form ~m ~n (h_rows_sparse : int array array) =
49
+
let k = n - m in
50
+
let nw = num_words n in
51
+
(* Build dense matrix from sparse rows *)
52
+
let mat = Array.init m (fun _ -> Array.make nw 0) in
53
+
Array.iteri
54
+
(fun i cols -> Array.iter (fun j -> mat_set mat i j 1) cols)
55
+
h_rows_sparse;
56
+
(* Column permutation: perm.(i) = original column at virtual position i.
57
+
We avoid physical column swaps by working with a virtual column mapping.
58
+
However, for the actual bit operations we need physical layout, so we
59
+
do swap the bits but in bulk (word at a time where possible). *)
60
+
let perm = Array.init n (fun i -> i) in
61
+
(* Swap two columns in the dense matrix and permutation *)
62
+
let swap_cols c1 c2 =
63
+
if c1 <> c2 then begin
64
+
let tmp = perm.(c1) in
65
+
perm.(c1) <- perm.(c2);
66
+
perm.(c2) <- tmp;
67
+
let w1 = c1 / word_bits and b1 = c1 mod word_bits in
68
+
let w2 = c2 / word_bits and b2 = c2 mod word_bits in
69
+
for r = 0 to m - 1 do
70
+
let row = mat.(r) in
71
+
let v1 = (row.(w1) lsr b1) land 1 in
72
+
let v2 = (row.(w2) lsr b2) land 1 in
73
+
if v1 <> v2 then begin
74
+
row.(w1) <- row.(w1) lxor (1 lsl b1);
75
+
row.(w2) <- row.(w2) lxor (1 lsl b2)
76
+
end
77
+
done
78
+
end
79
+
in
80
+
(* Put identity in the last m columns (columns k..n-1) by pivoting *)
81
+
for pivot = 0 to m - 1 do
82
+
let col = k + pivot in
83
+
(* Find a row >= pivot with a 1 in column col *)
84
+
let found_row = ref (-1) in
85
+
let r = ref pivot in
86
+
while !r < m && !found_row < 0 do
87
+
if mat_get mat !r col = 1 then found_row := !r;
88
+
incr r
89
+
done;
90
+
if !found_row < 0 then begin
91
+
(* No pivot in target column -- find a column in 0..k-1 that has
92
+
a 1 in some row >= pivot and swap it with col *)
93
+
let c = ref 0 in
94
+
while !c < k && !found_row < 0 do
95
+
r := pivot;
96
+
while !r < m && !found_row < 0 do
97
+
if mat_get mat !r !c = 1 then begin
98
+
swap_cols !c col;
99
+
found_row := !r
100
+
end;
101
+
incr r
102
+
done;
103
+
incr c
104
+
done
105
+
end;
106
+
if !found_row >= 0 then begin
107
+
(* Swap rows *)
108
+
if !found_row <> pivot then begin
109
+
let tmp = mat.(pivot) in
110
+
mat.(pivot) <- mat.(!found_row);
111
+
mat.(!found_row) <- tmp
112
+
end;
113
+
(* Eliminate: XOR pivot row into all other rows that have a 1 in col *)
114
+
for r2 = 0 to m - 1 do
115
+
if r2 <> pivot && mat_get mat r2 col = 1 then
116
+
mat_xor_row mat r2 pivot nw
117
+
done
118
+
end
119
+
done;
120
+
(* Extract generator parity sub-matrix:
121
+
gen_p.(j) = { i : 0 <= i < k, mat[j][i] = 1 } *)
122
+
let gen_p =
123
+
Array.init m (fun j ->
124
+
let cols = ref [] in
125
+
for i = k - 1 downto 0 do
126
+
if mat_get mat j i = 1 then cols := i :: !cols
127
+
done;
128
+
Array.of_list !cols)
129
+
in
130
+
(perm, gen_p)
131
+
132
+
(* --- Pseudo-random regular LDPC matrix construction --- *)
133
+
134
+
let make_regular_h ~m ~n ~wc ~wr:_ seed =
135
+
let state = ref seed in
136
+
let next_int bound =
137
+
state := ((!state * 1103515245) + 12345) land 0x3FFFFFFF;
138
+
!state mod bound
139
+
in
140
+
let h = Array.init m (fun _ -> Hashtbl.create 8) in
141
+
let add_edge r c =
142
+
if Hashtbl.mem h.(r) c then false
143
+
else begin
144
+
Hashtbl.replace h.(r) c true;
145
+
true
146
+
end
147
+
in
148
+
for sub = 0 to wc - 1 do
149
+
if sub = 0 then begin
150
+
for j = 0 to n - 1 do
151
+
ignore (add_edge (j mod m) j)
152
+
done
153
+
end
154
+
else begin
155
+
let perm = Array.init n (fun i -> i) in
156
+
for i = n - 1 downto 1 do
157
+
let j = next_int (i + 1) in
158
+
let tmp = perm.(i) in
159
+
perm.(i) <- perm.(j);
160
+
perm.(j) <- tmp
161
+
done;
162
+
for j = 0 to n - 1 do
163
+
let r = perm.(j) mod m in
164
+
let added = ref false in
165
+
for dr = 0 to m - 1 do
166
+
if not !added then begin
167
+
let r' = (r + dr) mod m in
168
+
if add_edge r' j then added := true
169
+
end
170
+
done
171
+
done
172
+
end
173
+
done;
174
+
let h_rows =
175
+
Array.init m (fun i ->
176
+
let cols = Hashtbl.fold (fun c _ acc -> c :: acc) h.(i) [] in
177
+
Array.of_list (List.sort compare cols))
178
+
in
179
+
let col_lists = Array.make n [] in
180
+
Array.iteri
181
+
(fun r cols ->
182
+
Array.iter (fun c -> col_lists.(c) <- r :: col_lists.(c)) cols)
183
+
h_rows;
184
+
let h_cols =
185
+
Array.map (fun rows -> Array.of_list (List.sort compare rows)) col_lists
186
+
in
187
+
(h_rows, h_cols)
188
+
189
+
let make_code ~n ~k ~m h_rows _h_cols =
190
+
let perm, gen_p = systematic_form ~m ~n h_rows in
191
+
(* Rebuild h_rows and h_cols with permuted columns *)
192
+
let inv_perm = Array.make n 0 in
193
+
Array.iteri (fun new_pos orig -> inv_perm.(orig) <- new_pos) perm;
194
+
let h_rows' =
195
+
Array.map
196
+
(fun cols ->
197
+
let cols' = Array.map (fun c -> inv_perm.(c)) cols in
198
+
Array.sort compare cols';
199
+
cols')
200
+
h_rows
201
+
in
202
+
let col_lists = Array.make n [] in
203
+
Array.iteri
204
+
(fun r cols ->
205
+
Array.iter (fun c -> col_lists.(c) <- r :: col_lists.(c)) cols)
206
+
h_rows';
207
+
let h_cols' =
208
+
Array.map (fun rows -> Array.of_list (List.sort compare rows)) col_lists
209
+
in
210
+
{ n; k; m; h_rows = h_rows'; h_cols = h_cols'; gen_p }
211
+
212
+
(* --- Rate 1/2 k=1024 code (lazy to avoid startup cost) --- *)
213
+
214
+
let ccsds_rate_1_2 =
215
+
let code =
216
+
lazy begin
217
+
let n = 2048 in
218
+
let k = 1024 in
219
+
let m = n - k in
220
+
let wc = 3 in
221
+
let wr = 6 in
222
+
let h_rows, h_cols = make_regular_h ~m ~n ~wc ~wr 42 in
223
+
make_code ~n ~k ~m h_rows h_cols
224
+
end
225
+
in
226
+
Lazy.force code
227
+
228
+
(* --- Bit manipulation helpers --- *)
229
+
230
+
let get_bit (buf : bytes) i =
231
+
let byte_idx = i / 8 in
232
+
let bit_pos = 7 - (i mod 8) in
233
+
(Char.code (Bytes.get buf byte_idx) lsr bit_pos) land 1
234
+
235
+
let set_bit (buf : bytes) i v =
236
+
let byte_idx = i / 8 in
237
+
let bit_pos = 7 - (i mod 8) in
238
+
let old = Char.code (Bytes.get buf byte_idx) in
239
+
let cleared = old land lnot (1 lsl bit_pos) in
240
+
Bytes.set buf byte_idx (Char.chr (cleared lor ((v land 1) lsl bit_pos)))
241
+
242
+
(* --- Systematic encoder --- *)
243
+
244
+
let encode code data =
245
+
let data_bits = Bytes.length data * 8 in
246
+
if data_bits < code.k then
247
+
invalid_arg
248
+
(Printf.sprintf
249
+
"Ldpc.encode: data too short, need %d bits (%d bytes), got %d bits"
250
+
code.k
251
+
((code.k + 7) / 8)
252
+
data_bits);
253
+
(* Read k information bits from data *)
254
+
let info = Array.init code.k (fun i -> get_bit data i) in
255
+
(* Compute m parity bits: parity_j = XOR of info_i for i in gen_p.(j) *)
256
+
let parity =
257
+
Array.init code.m (fun j ->
258
+
Array.fold_left (fun acc i -> acc lxor info.(i)) 0 code.gen_p.(j))
259
+
in
260
+
(* Build codeword: [info | parity], packed into bytes *)
261
+
let cw_bytes = (code.n + 7) / 8 in
262
+
let cw = Bytes.make cw_bytes '\x00' in
263
+
for i = 0 to code.k - 1 do
264
+
set_bit cw i info.(i)
265
+
done;
266
+
for j = 0 to code.m - 1 do
267
+
set_bit cw (code.k + j) parity.(j)
268
+
done;
269
+
cw
270
+
271
+
(* --- Syndrome check --- *)
272
+
273
+
let syndrome code (bits : int array) =
274
+
(* Returns true if all parity checks are satisfied *)
275
+
let ok = ref true in
276
+
for i = 0 to code.m - 1 do
277
+
let s =
278
+
Array.fold_left (fun acc j -> acc lxor bits.(j)) 0 code.h_rows.(i)
279
+
in
280
+
if s <> 0 then ok := false
281
+
done;
282
+
!ok
283
+
284
+
(* --- Sum-product (belief propagation) decoder --- *)
285
+
286
+
(* LLR sign convention: positive LLR means bit is more likely 0,
287
+
negative means more likely 1.
288
+
289
+
For hard-decision input we model a BSC with crossover probability ~0.045,
290
+
giving channel LLR magnitude = ln((1-p)/p) ~ 3.0. This provides enough
291
+
soft information for BP to converge while remaining realistic for the
292
+
error rates LDPC codes are designed to handle. *)
293
+
let hard_decision_llr = 3.0
294
+
295
+
let decode ?(max_iter = 50) code codeword =
296
+
let cw_bits = Bytes.length codeword * 8 in
297
+
if cw_bits < code.n then
298
+
Error
299
+
(Printf.sprintf
300
+
"Ldpc.decode: codeword too short, need %d bits (%d bytes), got %d bits"
301
+
code.n
302
+
((code.n + 7) / 8)
303
+
cw_bits)
304
+
else begin
305
+
(* Initialize LLRs from hard-decision input *)
306
+
let channel_llr =
307
+
Array.init code.n (fun i ->
308
+
if get_bit codeword i = 0 then hard_decision_llr
309
+
else Float.neg hard_decision_llr)
310
+
in
311
+
(* Message arrays.
312
+
For each edge (check i, variable j):
313
+
- r.(i).(local_j) = check-to-variable message
314
+
- q.(i).(local_j) = variable-to-check message
315
+
316
+
We index edges by (check_node, position_in_check_node_neighbors). *)
317
+
let num_checks = code.m in
318
+
319
+
(* For each check node i, h_rows.(i) gives the connected variable nodes.
320
+
For each variable node j, h_cols.(j) gives the connected check nodes. *)
321
+
322
+
(* check-to-variable messages, indexed [check][local_var_idx] *)
323
+
let r =
324
+
Array.init num_checks (fun i ->
325
+
Array.make (Array.length code.h_rows.(i)) 0.0)
326
+
in
327
+
(* variable-to-check messages, indexed [check][local_var_idx] *)
328
+
let q =
329
+
Array.init num_checks (fun i ->
330
+
Array.init
331
+
(Array.length code.h_rows.(i))
332
+
(fun local_j ->
333
+
let j = code.h_rows.(i).(local_j) in
334
+
channel_llr.(j)))
335
+
in
336
+
(* Build reverse index: for variable j and check i, what is the local
337
+
index of j in h_rows.(i)? *)
338
+
let var_to_check_local =
339
+
Array.init code.n (fun j ->
340
+
Array.map
341
+
(fun i ->
342
+
(* Find local index of j in h_rows.(i) *)
343
+
let local = ref 0 in
344
+
let found = ref false in
345
+
for idx = 0 to Array.length code.h_rows.(i) - 1 do
346
+
if (not !found) && code.h_rows.(i).(idx) = j then begin
347
+
local := idx;
348
+
found := true
349
+
end
350
+
done;
351
+
(i, !local))
352
+
code.h_cols.(j))
353
+
in
354
+
355
+
(* Posterior LLR for each variable node *)
356
+
let posterior = Array.copy channel_llr in
357
+
358
+
let converged = ref false in
359
+
let iter = ref 0 in
360
+
361
+
while (not !converged) && !iter < max_iter do
362
+
incr iter;
363
+
364
+
(* --- Check node update (horizontal step) --- *)
365
+
for i = 0 to num_checks - 1 do
366
+
let deg = Array.length code.h_rows.(i) in
367
+
let signs =
368
+
Array.init deg (fun lj -> if q.(i).(lj) >= 0.0 then 1 else -1)
369
+
in
370
+
let phi x =
371
+
let ax = Float.abs x in
372
+
if ax < 1e-10 then 20.0 (* clamp *)
373
+
else if ax > 19.0 then 1e-9
374
+
else
375
+
let e = exp ax in
376
+
Float.abs (log ((e +. 1.0) /. (e -. 1.0)))
377
+
in
378
+
let phis = Array.init deg (fun lj -> phi q.(i).(lj)) in
379
+
let total_sign = Array.fold_left ( * ) 1 signs in
380
+
let total_phi = Array.fold_left ( +. ) 0.0 phis in
381
+
for lj = 0 to deg - 1 do
382
+
let excl_sign = total_sign * signs.(lj) in
383
+
let excl_phi = total_phi -. phis.(lj) in
384
+
let excl_phi = Float.max excl_phi 1e-10 in
385
+
let magnitude = phi excl_phi in
386
+
r.(i).(lj) <- (if excl_sign > 0 then 1.0 else -1.0) *. magnitude
387
+
done
388
+
done;
389
+
390
+
(* --- Variable node update (vertical step) --- *)
391
+
for j = 0 to code.n - 1 do
392
+
let checks = var_to_check_local.(j) in
393
+
let num_checks_j = Array.length checks in
394
+
let total_r = ref channel_llr.(j) in
395
+
for idx = 0 to num_checks_j - 1 do
396
+
let i, lj = checks.(idx) in
397
+
total_r := !total_r +. r.(i).(lj)
398
+
done;
399
+
posterior.(j) <- !total_r;
400
+
for idx = 0 to num_checks_j - 1 do
401
+
let i, lj = checks.(idx) in
402
+
q.(i).(lj) <- !total_r -. r.(i).(lj)
403
+
done
404
+
done;
405
+
406
+
(* --- Hard decision and syndrome check --- *)
407
+
let hard =
408
+
Array.init code.n (fun j -> if posterior.(j) >= 0.0 then 0 else 1)
409
+
in
410
+
if syndrome code hard then converged := true
411
+
done;
412
+
413
+
(* Extract the decoded information bits *)
414
+
let hard =
415
+
Array.init code.n (fun j -> if posterior.(j) >= 0.0 then 0 else 1)
416
+
in
417
+
let data_bytes = (code.k + 7) / 8 in
418
+
let out = Bytes.make data_bytes '\x00' in
419
+
for i = 0 to code.k - 1 do
420
+
set_bit out i hard.(i)
421
+
done;
422
+
Ok out
423
+
end
+41
lib/ldpc.mli
+41
lib/ldpc.mli
···
1
+
(** Low-Density Parity-Check codes with belief propagation decoding.
2
+
3
+
Generic LDPC encoder/decoder framework with systematic encoding and
4
+
sum-product (belief propagation) decoding. Configurable parity check matrix.
5
+
6
+
@see <https://public.ccsds.org/Pubs/131x0b4.pdf> CCSDS 131.0-B-4 Section 7
7
+
@see <https://public.ccsds.org/Pubs/231x0b4.pdf> CCSDS 231.0-B Section 4 *)
8
+
9
+
type code
10
+
(** An LDPC code defined by its parity check matrix, generator matrix, and code
11
+
parameters (n, k, m). *)
12
+
13
+
val ccsds_rate_1_2 : code
14
+
(** Rate 1/2 LDPC code with k=1024, n=2048. Uses a pseudo-random regular (3,6)
15
+
parity check matrix for testing the encoder/decoder framework.
16
+
17
+
{b NOT spec-compliant:} this does not use the CCSDS 131.0-B-4 Annex G
18
+
parity-check matrix. The actual CCSDS codes use circulant-based H matrices
19
+
defined by base matrices in Annex G Tables G-1 through G-8. Codewords
20
+
produced by this code are not interoperable with CCSDS-compliant decoders.
21
+
*)
22
+
23
+
val encode : code -> bytes -> bytes
24
+
(** [encode code data] performs systematic LDPC encoding.
25
+
26
+
The input [data] must contain at least [k] bits (i.e., [ceil(k/8)] bytes).
27
+
The output is a codeword of [ceil(n/8)] bytes containing the [k] information
28
+
bits followed by [n-k] parity bits, packed MSB-first.
29
+
30
+
@raise Invalid_argument if the input is too short. *)
31
+
32
+
val decode : ?max_iter:int -> code -> bytes -> (bytes, string) result
33
+
(** [decode ?max_iter code codeword] decodes an LDPC codeword using the
34
+
sum-product (belief propagation) algorithm.
35
+
36
+
The input [codeword] must contain at least [n] bits. Hard-decision input is
37
+
converted to LLRs (bit 0 -> +1.0, bit 1 -> -1.0). The decoder iterates until
38
+
the syndrome is zero or [max_iter] iterations are reached (default 50).
39
+
40
+
Returns [Ok data] with the [ceil(k/8)] decoded information bytes on success,
41
+
or [Error msg] if the input is too short. *)
+1
test/test.ml
+1
test/test.ml
···
1
+
let () = Alcotest.run "ldpc" [ Test_ldpc.suite ]
+116
test/test_ldpc.ml
+116
test/test_ldpc.ml
···
1
+
(** Tests for LDPC codes with belief propagation decoding.
2
+
3
+
- {b Roundtrip}: encode then decode with no channel errors.
4
+
- {b Syndrome}: verify that a valid codeword satisfies H*c = 0.
5
+
- {b Error correction}: flip 5 bits in a codeword, verify BP decoding
6
+
recovers the original data.
7
+
- {b All-zeros / all-ones}: edge cases for systematic encoding.
8
+
- {b Codeword length}: verify n = 2048 bits (256 bytes) for rate 1/2.
9
+
- {b Short input}: verify decode rejects undersized codewords. *)
10
+
11
+
let ldpc = Ldpc.ccsds_rate_1_2
12
+
13
+
(** Helper: flip bit i in a bytes buffer. *)
14
+
let flip_bit buf i =
15
+
let byte_idx = i / 8 in
16
+
let bit_pos = 7 - (i mod 8) in
17
+
let old = Char.code (Bytes.get buf byte_idx) in
18
+
Bytes.set buf byte_idx (Char.chr (old lxor (1 lsl bit_pos)))
19
+
20
+
(** LDPC encode/decode roundtrip with no errors. *)
21
+
let test_ldpc_roundtrip () =
22
+
let data = Bytes.init 128 (fun i -> Char.chr (((i * 37) + 13) land 0xFF)) in
23
+
let codeword = Ldpc.encode ldpc data in
24
+
Alcotest.(check int) "LDPC codeword length" 256 (Bytes.length codeword);
25
+
Alcotest.(check string)
26
+
"LDPC systematic: data preserved" (Bytes.to_string data)
27
+
(Bytes.sub_string codeword 0 128);
28
+
match Ldpc.decode ldpc codeword with
29
+
| Ok recovered ->
30
+
Alcotest.(check string)
31
+
"LDPC roundtrip" (Bytes.to_string data)
32
+
(Bytes.to_string recovered)
33
+
| Error e -> Alcotest.fail (Printf.sprintf "LDPC decode failed: %s" e)
34
+
35
+
(** LDPC syndrome check: decoder converges immediately on an uncorrupted
36
+
codeword. *)
37
+
let test_ldpc_syndrome () =
38
+
let data = Bytes.init 128 (fun i -> Char.chr (((i * 53) + 7) land 0xFF)) in
39
+
let codeword = Ldpc.encode ldpc data in
40
+
match Ldpc.decode ~max_iter:1 ldpc codeword with
41
+
| Ok recovered ->
42
+
Alcotest.(check string)
43
+
"LDPC syndrome valid (1 iteration enough)" (Bytes.to_string data)
44
+
(Bytes.to_string recovered)
45
+
| Error e -> Alcotest.fail (Printf.sprintf "LDPC syndrome check failed: %s" e)
46
+
47
+
(** LDPC error correction: flip 5 bits and verify the decoder corrects them. *)
48
+
let test_ldpc_error_correction () =
49
+
let data = Bytes.init 128 (fun i -> Char.chr (((i * 41) + 99) land 0xFF)) in
50
+
let codeword = Ldpc.encode ldpc data in
51
+
let corrupted = Bytes.copy codeword in
52
+
let error_positions = [| 10; 200; 500; 1000; 1800 |] in
53
+
Array.iter (fun pos -> flip_bit corrupted pos) error_positions;
54
+
Alcotest.(check bool) "corrupted differs" true (corrupted <> codeword);
55
+
match Ldpc.decode ~max_iter:50 ldpc corrupted with
56
+
| Ok recovered ->
57
+
Alcotest.(check string)
58
+
"LDPC corrects 5 bit errors" (Bytes.to_string data)
59
+
(Bytes.to_string recovered)
60
+
| Error e ->
61
+
Alcotest.fail (Printf.sprintf "LDPC error correction failed: %s" e)
62
+
63
+
(** LDPC roundtrip with all-zeros data. *)
64
+
let test_ldpc_all_zeros () =
65
+
let data = Bytes.make 128 '\x00' in
66
+
let codeword = Ldpc.encode ldpc data in
67
+
let expected = Bytes.make 256 '\x00' in
68
+
Alcotest.(check string)
69
+
"LDPC all-zeros codeword" (Bytes.to_string expected)
70
+
(Bytes.to_string codeword);
71
+
match Ldpc.decode ldpc codeword with
72
+
| Ok recovered ->
73
+
Alcotest.(check string)
74
+
"LDPC all-zeros roundtrip" (Bytes.to_string data)
75
+
(Bytes.to_string recovered)
76
+
| Error e -> Alcotest.fail (Printf.sprintf "LDPC all-zeros failed: %s" e)
77
+
78
+
(** LDPC roundtrip with all-ones data. *)
79
+
let test_ldpc_all_ones () =
80
+
let data = Bytes.make 128 '\xFF' in
81
+
let codeword = Ldpc.encode ldpc data in
82
+
Alcotest.(check int)
83
+
"LDPC all-ones codeword length" 256 (Bytes.length codeword);
84
+
match Ldpc.decode ldpc codeword with
85
+
| Ok recovered ->
86
+
Alcotest.(check string)
87
+
"LDPC all-ones roundtrip" (Bytes.to_string data)
88
+
(Bytes.to_string recovered)
89
+
| Error e -> Alcotest.fail (Printf.sprintf "LDPC all-ones failed: %s" e)
90
+
91
+
(** LDPC codeword length check. *)
92
+
let test_ldpc_codeword_length () =
93
+
let data = Bytes.make 128 '\x42' in
94
+
let codeword = Ldpc.encode ldpc data in
95
+
Alcotest.(check int) "LDPC n=2048 -> 256 bytes" 256 (Bytes.length codeword)
96
+
97
+
(** LDPC decode with too-short input should return Error. *)
98
+
let test_ldpc_decode_short_input () =
99
+
let short = Bytes.make 100 '\x00' in
100
+
match Ldpc.decode ldpc short with
101
+
| Ok _ -> Alcotest.fail "LDPC decode should reject short input"
102
+
| Error _ -> ()
103
+
104
+
let suite =
105
+
( "ldpc",
106
+
[
107
+
Alcotest.test_case "LDPC roundtrip (no errors)" `Quick test_ldpc_roundtrip;
108
+
Alcotest.test_case "LDPC syndrome check" `Quick test_ldpc_syndrome;
109
+
Alcotest.test_case "LDPC error correction (5 bit flips)" `Quick
110
+
test_ldpc_error_correction;
111
+
Alcotest.test_case "LDPC all-zeros roundtrip" `Quick test_ldpc_all_zeros;
112
+
Alcotest.test_case "LDPC all-ones roundtrip" `Quick test_ldpc_all_ones;
113
+
Alcotest.test_case "LDPC codeword length" `Quick test_ldpc_codeword_length;
114
+
Alcotest.test_case "LDPC decode short input" `Quick
115
+
test_ldpc_decode_short_input;
116
+
] )