CCSDS 131.4-B short block-length LDPC codes
short-ldpc: extract decode/systematic_form helpers, add fuzz mli
- Split decode (104 lines) into phi, check_node_update, var_node_update,
build_var_to_check_local helpers; extract err_codeword_too_short.
- Split systematic_form into pivot_row, pivot_via_swap, eliminate_pivot,
swap_cols_fn helpers.
- Split regular_h into add_permuted_layer and extract_sparse_h helpers.
- Rename fuzz suite "short-ldpc" -> "short_ldpc" to match filename.
- Add fuzz_short_ldpc.mli.
+192
-148
+1
-1
fuzz/fuzz_short_ldpc.ml
+1
-1
fuzz/fuzz_short_ldpc.ml
+5
fuzz/fuzz_short_ldpc.mli
+5
fuzz/fuzz_short_ldpc.mli
+186
-147
lib/short_ldpc.ml
+186
-147
lib/short_ldpc.ml
···
40
40
d.(w) <- d.(w) lxor s.(w)
41
41
done
42
42
43
+
(* Locate a row >= [pivot] with a 1 in column [col]; returns -1 if none. *)
44
+
let pivot_row mat ~pivot ~col ~m =
45
+
let r = ref pivot in
46
+
let found = ref (-1) in
47
+
while !r < m && !found < 0 do
48
+
if mat_get mat !r col = 1 then found := !r;
49
+
incr r
50
+
done;
51
+
!found
52
+
53
+
(* Search columns 0..k-1 for a row >= [pivot] with a 1, swap that column
54
+
into [col] via [swap_cols], and return the found row index (or -1). *)
55
+
let pivot_via_swap mat ~swap_cols ~pivot ~col ~k ~m =
56
+
let c = ref 0 in
57
+
let found = ref (-1) in
58
+
while !c < k && !found < 0 do
59
+
let r = ref pivot in
60
+
while !r < m && !found < 0 do
61
+
if mat_get mat !r !c = 1 then begin
62
+
swap_cols !c col;
63
+
found := !r
64
+
end;
65
+
incr r
66
+
done;
67
+
incr c
68
+
done;
69
+
!found
70
+
43
71
(* Gaussian elimination over GF(2) to put H into [A | I_m] form. *)
72
+
(* Build a column-swap closure: swaps columns c1/c2 in [perm] and [mat]. *)
73
+
let swap_cols_fn ~m perm mat =
74
+
fun c1 c2 ->
75
+
if c1 <> c2 then begin
76
+
let tmp = perm.(c1) in
77
+
perm.(c1) <- perm.(c2);
78
+
perm.(c2) <- tmp;
79
+
let w1 = c1 / word_bits and b1 = c1 mod word_bits in
80
+
let w2 = c2 / word_bits and b2 = c2 mod word_bits in
81
+
for r = 0 to m - 1 do
82
+
let row = mat.(r) in
83
+
let v1 = (row.(w1) lsr b1) land 1 in
84
+
let v2 = (row.(w2) lsr b2) land 1 in
85
+
if v1 <> v2 then begin
86
+
row.(w1) <- row.(w1) lxor (1 lsl b1);
87
+
row.(w2) <- row.(w2) lxor (1 lsl b2)
88
+
end
89
+
done
90
+
end
91
+
92
+
(* Eliminate one pivot column: bring [found_row] to position [pivot] then
93
+
XOR it into every other row that has a 1 in [col]. *)
94
+
let eliminate_pivot ~m ~nw mat ~pivot ~col ~found_row =
95
+
if found_row <> pivot then begin
96
+
let tmp = mat.(pivot) in
97
+
mat.(pivot) <- mat.(found_row);
98
+
mat.(found_row) <- tmp
99
+
end;
100
+
for r2 = 0 to m - 1 do
101
+
if r2 <> pivot && mat_get mat r2 col = 1 then mat_xor_row mat r2 pivot nw
102
+
done
103
+
44
104
let systematic_form ~m ~n (h_rows_sparse : int array array) =
45
105
let k = n - m in
46
106
let nw = num_words n in
···
49
109
(fun i cols -> Array.iter (fun j -> mat_set mat i j 1) cols)
50
110
h_rows_sparse;
51
111
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
112
+
let swap_cols = swap_cols_fn ~m perm mat in
70
113
for pivot = 0 to m - 1 do
71
114
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
115
+
let found_row =
116
+
let r = pivot_row mat ~pivot ~col ~m in
117
+
if r >= 0 then r else pivot_via_swap mat ~swap_cols ~pivot ~col ~k ~m
118
+
in
119
+
if found_row >= 0 then eliminate_pivot ~m ~nw mat ~pivot ~col ~found_row
103
120
done;
104
121
let gen_p =
105
122
Array.init m (fun j ->
···
113
130
114
131
(* --- Pseudo-random regular LDPC matrix construction --- *)
115
132
133
+
(* Add a permuted layer of edges to [h]: each variable picks a row from a
134
+
shuffled permutation, scanning forward on collision until an edge can be
135
+
added. *)
136
+
let add_permuted_layer ~m ~n ~next_int ~add_edge =
137
+
let perm = Array.init n (fun i -> i) in
138
+
for i = n - 1 downto 1 do
139
+
let j = next_int (i + 1) in
140
+
let tmp = perm.(i) in
141
+
perm.(i) <- perm.(j);
142
+
perm.(j) <- tmp
143
+
done;
144
+
for j = 0 to n - 1 do
145
+
let r = perm.(j) mod m in
146
+
let added = ref false in
147
+
for dr = 0 to m - 1 do
148
+
if not !added then begin
149
+
let r' = (r + dr) mod m in
150
+
if add_edge r' j then added := true
151
+
end
152
+
done
153
+
done
154
+
155
+
(* Extract sorted row/column adjacency arrays from a hashtable-backed sparse
156
+
matrix [h]. *)
157
+
let extract_sparse_h ~m ~n h =
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
+
116
173
(* For short codes we use lighter column weights to keep the graph sparse
117
174
enough for BP to converge on small blocks. *)
118
175
let regular_h ~m ~n ~wc ~wr:_ seed =
···
130
187
end
131
188
in
132
189
for sub = 0 to wc - 1 do
133
-
if sub = 0 then begin
190
+
if sub = 0 then
134
191
for j = 0 to n - 1 do
135
192
ignore (add_edge (j mod m) j)
136
193
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
194
+
else add_permuted_layer ~m ~n ~next_int ~add_edge
157
195
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)
196
+
extract_sparse_h ~m ~n h
172
197
173
198
let code ~n ~k ~m h_rows _h_cols =
174
199
let perm, gen_p = systematic_form ~m ~n h_rows in
···
293
318
294
319
let hard_decision_llr = 3.0
295
320
321
+
let err_codeword_too_short ~n ~cw_bits =
322
+
Error
323
+
(Fmt.str
324
+
"Short_ldpc.decode: codeword too short, need %d bits (%d bytes), got %d \
325
+
bits"
326
+
n
327
+
((n + 7) / 8)
328
+
cw_bits)
329
+
330
+
(* Phi function used in the sum-product check-node update. *)
331
+
let phi x =
332
+
let ax = Float.abs x in
333
+
if ax < 1e-10 then 20.0
334
+
else if ax > 19.0 then 1e-9
335
+
else
336
+
let e = exp ax in
337
+
Float.abs (log ((e +. 1.0) /. (e -. 1.0)))
338
+
339
+
(* Check-node message update for one iteration: r.(i).(lj) is the message
340
+
from check i to variable [code.h_rows.(i).(lj)] excluding that variable. *)
341
+
let check_node_update ~num_checks ~code ~q ~r =
342
+
for i = 0 to num_checks - 1 do
343
+
let deg = Array.length code.h_rows.(i) in
344
+
let signs =
345
+
Array.init deg (fun lj -> if q.(i).(lj) >= 0.0 then 1 else -1)
346
+
in
347
+
let phis = Array.init deg (fun lj -> phi q.(i).(lj)) in
348
+
let total_sign = Array.fold_left ( * ) 1 signs in
349
+
let total_phi = Array.fold_left ( +. ) 0.0 phis in
350
+
for lj = 0 to deg - 1 do
351
+
let excl_sign = total_sign * signs.(lj) in
352
+
let excl_phi = total_phi -. phis.(lj) in
353
+
let excl_phi = Float.max excl_phi 1e-10 in
354
+
let magnitude = phi excl_phi in
355
+
r.(i).(lj) <- (if excl_sign > 0 then 1.0 else -1.0) *. magnitude
356
+
done
357
+
done
358
+
359
+
(* Variable-node message update: aggregates r messages back into q messages
360
+
and writes the per-variable posterior LLR. *)
361
+
let var_node_update ~code ~channel_llr ~var_to_check_local ~r ~q ~posterior =
362
+
for j = 0 to code.n - 1 do
363
+
let checks = var_to_check_local.(j) in
364
+
let num_checks_j = Array.length checks in
365
+
let total_r = ref channel_llr.(j) in
366
+
for idx = 0 to num_checks_j - 1 do
367
+
let i, lj = checks.(idx) in
368
+
total_r := !total_r +. r.(i).(lj)
369
+
done;
370
+
posterior.(j) <- !total_r;
371
+
for idx = 0 to num_checks_j - 1 do
372
+
let i, lj = checks.(idx) in
373
+
q.(i).(lj) <- !total_r -. r.(i).(lj)
374
+
done
375
+
done
376
+
377
+
(* Reverse mapping: for each variable j, list (check_i, local_index) pairs
378
+
identifying where j appears in code.h_rows.(check_i). *)
379
+
let build_var_to_check_local code =
380
+
Array.init code.n (fun j ->
381
+
Array.map
382
+
(fun i ->
383
+
let local = ref 0 in
384
+
let found = ref false in
385
+
for idx = 0 to Array.length code.h_rows.(i) - 1 do
386
+
if (not !found) && code.h_rows.(i).(idx) = j then begin
387
+
local := idx;
388
+
found := true
389
+
end
390
+
done;
391
+
(i, !local))
392
+
code.h_cols.(j))
393
+
296
394
let decode ?(max_iter = 50) code codeword =
297
395
let cw_bits = Bytes.length codeword * 8 in
298
-
if cw_bits < code.n then
299
-
Error
300
-
(Fmt.str
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)
396
+
if cw_bits < code.n then err_codeword_too_short ~n:code.n ~cw_bits
306
397
else begin
307
398
let channel_llr =
308
399
Array.init code.n (fun i ->
···
322
413
let j = code.h_rows.(i).(local_j) in
323
414
channel_llr.(j)))
324
415
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
416
+
let var_to_check_local = build_var_to_check_local code in
340
417
let posterior = Array.copy channel_llr in
341
418
let converged = ref false in
342
419
let iter = ref 0 in
343
420
while (not !converged) && !iter < max_iter do
344
421
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;
422
+
check_node_update ~num_checks ~code ~q ~r;
423
+
var_node_update ~code ~channel_llr ~var_to_check_local ~r ~q ~posterior;
385
424
let hard =
386
425
Array.init code.n (fun j -> if posterior.(j) >= 0.0 then 0 else 1)
387
426
in