LDPC codes with belief propagation decoding
Add --here, --location CITY, and Arg.enum to contact CLI
Location can be specified three ways:
--lat 34.05 --lon=-118.25 (coordinates)
--location la (city preset via Arg.enum)
--here (auto-detect via IP geolocation)
Known cities: la, sf, nyc, london, paris, tokyo, ksc, kourou,
baikonur, bangalore, sydney. Cmdliner provides tab completion
and error messages for invalid city names.
+426
-86
+410
-75
lib/ldpc.ml
+410
-75
lib/ldpc.ml
···
1
1
(** Low-Density Parity-Check codes with belief propagation decoding.
2
2
3
-
Generic LDPC encoder/decoder framework with systematic encoding and
4
-
sum-product (belief propagation) decoding. *)
3
+
Implements the CCSDS 131.0-B-4 AR4JA (Accumulate-Repeat-4-Jagged-
4
+
Accumulate) codes with circulant-based parity-check matrices, systematic
5
+
encoding, and sum-product belief propagation decoding. *)
5
6
6
7
(* Sparse representation of the parity check matrix H.
7
8
Each row is a list of column indices where H has a 1.
8
9
Each column is a list of row indices where H has a 1. *)
9
10
type code = {
10
-
n : int; (* codeword length *)
11
+
n : int; (* transmitted codeword length *)
11
12
k : int; (* information length *)
12
-
m : int; (* number of parity check equations = n - k *)
13
+
m : int; (* number of parity check equations *)
14
+
n_full : int; (* full codeword length including punctured bits *)
15
+
punctured : int; (* number of punctured bits *)
13
16
h_rows : int array array; (* h_rows.(i) = columns with 1 in row i *)
14
17
h_cols : int array array; (* h_cols.(j) = rows with 1 in column j *)
15
18
gen_p :
···
129
132
in
130
133
(perm, gen_p)
131
134
132
-
(* --- Pseudo-random regular LDPC matrix construction --- *)
135
+
(* --- CCSDS AR4JA parity-check matrix construction ---
136
+
137
+
The AR4JA codes use a protograph-based construction with a two-step
138
+
circulant expansion. The protograph for rate 1/2 has 3 check nodes,
139
+
5 variable nodes (4 transmitted + 1 punctured), giving effective rate
140
+
k/(n_transmitted) = 1/2.
141
+
142
+
The first expansion by factor 4 produces a 4x11 intermediate weight
143
+
matrix. The second expansion by factor M/4 (circulant_size) produces
144
+
the full H matrix.
145
+
146
+
Data from CCSDS 131.0-B-4 / 131.1-O-2, verified against the
147
+
labrador-ldpc reference implementation.
148
+
149
+
Reference: CCSDS 131.0-B-4, TM Synchronization and Channel Coding *)
150
+
151
+
(* Block types in the prototype matrix *)
152
+
let _hz = 0 (* zero block *)
153
+
let hi = 1 (* identity block (possibly shifted) *)
154
+
let hp = 2 (* permutation block using theta/phi tables *)
155
+
156
+
(* Prototype matrix for rate 1/2: 3 used rows x 11 columns.
157
+
Each entry is (block_type, shift_value).
158
+
The 4th row is all zeros (unused).
159
+
Three layers are superimposed (up to 3 edges per protograph edge). *)
160
+
161
+
(* Layer 0 *)
162
+
let tm_r12_proto_0 =
163
+
[|
164
+
(* row 0 *)
165
+
[|
166
+
(_hz, 0);
167
+
(_hz, 0);
168
+
(hi, 0);
169
+
(_hz, 0);
170
+
(hi, 0);
171
+
(0, 0);
172
+
(0, 0);
173
+
(0, 0);
174
+
(0, 0);
175
+
(0, 0);
176
+
(0, 0);
177
+
|];
178
+
(* row 1 *)
179
+
[|
180
+
(hi, 0);
181
+
(hi, 0);
182
+
(_hz, 0);
183
+
(hi, 0);
184
+
(hp, 1);
185
+
(0, 0);
186
+
(0, 0);
187
+
(0, 0);
188
+
(0, 0);
189
+
(0, 0);
190
+
(0, 0);
191
+
|];
192
+
(* row 2 *)
193
+
[|
194
+
(hi, 0);
195
+
(hp, 4);
196
+
(_hz, 0);
197
+
(hp, 6);
198
+
(hi, 0);
199
+
(0, 0);
200
+
(0, 0);
201
+
(0, 0);
202
+
(0, 0);
203
+
(0, 0);
204
+
(0, 0);
205
+
|];
206
+
|]
207
+
208
+
(* Layer 1 *)
209
+
let tm_r12_proto_1 =
210
+
[|
211
+
(* row 0 *)
212
+
[|
213
+
(0, 0);
214
+
(0, 0);
215
+
(0, 0);
216
+
(0, 0);
217
+
(hp, 0);
218
+
(0, 0);
219
+
(0, 0);
220
+
(0, 0);
221
+
(0, 0);
222
+
(0, 0);
223
+
(0, 0);
224
+
|];
225
+
(* row 1 *)
226
+
[|
227
+
(0, 0);
228
+
(0, 0);
229
+
(0, 0);
230
+
(0, 0);
231
+
(hp, 2);
232
+
(0, 0);
233
+
(0, 0);
234
+
(0, 0);
235
+
(0, 0);
236
+
(0, 0);
237
+
(0, 0);
238
+
|];
239
+
(* row 2 *)
240
+
[|
241
+
(0, 0);
242
+
(hp, 5);
243
+
(0, 0);
244
+
(hp, 7);
245
+
(0, 0);
246
+
(0, 0);
247
+
(0, 0);
248
+
(0, 0);
249
+
(0, 0);
250
+
(0, 0);
251
+
(0, 0);
252
+
|];
253
+
|]
254
+
255
+
(* Layer 2 *)
256
+
let tm_r12_proto_2 =
257
+
[|
258
+
(* row 0 *)
259
+
[|
260
+
(0, 0);
261
+
(0, 0);
262
+
(0, 0);
263
+
(0, 0);
264
+
(0, 0);
265
+
(0, 0);
266
+
(0, 0);
267
+
(0, 0);
268
+
(0, 0);
269
+
(0, 0);
270
+
(0, 0);
271
+
|];
272
+
(* row 1 *)
273
+
[|
274
+
(0, 0);
275
+
(0, 0);
276
+
(0, 0);
277
+
(0, 0);
278
+
(hp, 3);
279
+
(0, 0);
280
+
(0, 0);
281
+
(0, 0);
282
+
(0, 0);
283
+
(0, 0);
284
+
(0, 0);
285
+
|];
286
+
(* row 2 *)
287
+
[|
288
+
(0, 0);
289
+
(0, 0);
290
+
(0, 0);
291
+
(0, 0);
292
+
(0, 0);
293
+
(0, 0);
294
+
(0, 0);
295
+
(0, 0);
296
+
(0, 0);
297
+
(0, 0);
298
+
(0, 0);
299
+
|];
300
+
|]
301
+
302
+
(* THETA_K: row permutation index, 26 entries (k=0..25) *)
303
+
let theta_k =
304
+
[|
305
+
3; 0; 1; 2; 2; 3; 0; 1; 0; 1; 2; 0; 2; 3; 0; 1; 2; 0; 1; 2; 0; 1; 2; 1; 2; 3;
306
+
|]
307
+
308
+
(* PHI_J_K for M=128 (circulant_size = M/4 = 32, so values are mod 32).
309
+
phi.(j).(k) for j in 0..3, k in 0..25.
310
+
Used for rate 1/2, k=1024 (submatrix_size M=512, circulant_size=128). *)
311
+
let phi_j_k_m128 =
312
+
[|
313
+
[|
314
+
1;
315
+
22;
316
+
0;
317
+
26;
318
+
0;
319
+
10;
320
+
5;
321
+
18;
322
+
3;
323
+
22;
324
+
3;
325
+
8;
326
+
25;
327
+
25;
328
+
2;
329
+
27;
330
+
7;
331
+
7;
332
+
15;
333
+
10;
334
+
4;
335
+
19;
336
+
7;
337
+
9;
338
+
26;
339
+
17;
340
+
|];
341
+
[|
342
+
0;
343
+
27;
344
+
30;
345
+
28;
346
+
7;
347
+
1;
348
+
8;
349
+
20;
350
+
26;
351
+
24;
352
+
4;
353
+
12;
354
+
23;
355
+
15;
356
+
15;
357
+
22;
358
+
31;
359
+
3;
360
+
29;
361
+
21;
362
+
2;
363
+
5;
364
+
11;
365
+
26;
366
+
9;
367
+
17;
368
+
|];
369
+
[|
370
+
0;
371
+
12;
372
+
30;
373
+
18;
374
+
10;
375
+
16;
376
+
13;
377
+
9;
378
+
7;
379
+
15;
380
+
16;
381
+
18;
382
+
4;
383
+
23;
384
+
5;
385
+
3;
386
+
29;
387
+
11;
388
+
4;
389
+
8;
390
+
2;
391
+
11;
392
+
11;
393
+
3;
394
+
15;
395
+
13;
396
+
|];
397
+
[|
398
+
0;
399
+
13;
400
+
19;
401
+
14;
402
+
15;
403
+
20;
404
+
17;
405
+
4;
406
+
4;
407
+
11;
408
+
17;
409
+
20;
410
+
8;
411
+
22;
412
+
19;
413
+
15;
414
+
5;
415
+
21;
416
+
17;
417
+
9;
418
+
20;
419
+
18;
420
+
31;
421
+
13;
422
+
2;
423
+
18;
424
+
|];
425
+
|]
426
+
427
+
(* Construct the CCSDS AR4JA parity-check matrix H for rate 1/2, k=1024.
133
428
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
429
+
The prototype is a 3x11 matrix (with 3 layers for superimposed edges).
430
+
Each block is M x M where M = submatrix_size = 512.
431
+
The resulting H matrix has dimensions:
432
+
rows = 3 * M = 1536
433
+
cols = 11 * M = 5632
434
+
But only the first 5*M = 2560 columns are non-trivial (columns 5..10
435
+
are all-zero in the prototype, so those blocks contribute nothing).
436
+
437
+
The full code has n_full = 2560, k = 1024, m = 1536.
438
+
After puncturing 512 bits (the last block column, index 4), the
439
+
transmitted codeword has n = 2048 bits.
440
+
441
+
For each non-zero block at prototype position (row_idx, col_idx):
442
+
- For each check index c in 0..M-1:
443
+
- HI block: H[row_idx*M + c, col_idx*M + ((c + shift) mod M)] = 1
444
+
- HP block: H[row_idx*M + c, col_idx*M + pi_k(c)] = 1
445
+
where pi_k uses THETA_K and PHI_J_K tables *)
446
+
447
+
let make_ar4ja_h () =
448
+
let m_sub = 512 in
449
+
(* submatrix size *)
450
+
let circ = m_sub / 4 in
451
+
(* circulant size = 128 *)
452
+
let log_circ = 7 in
453
+
(* log2(128) = 7 *)
454
+
let mod_m = m_sub - 1 in
455
+
(* bitmask for mod M *)
456
+
let mod_circ = circ - 1 in
457
+
(* bitmask for mod circulant_size *)
458
+
let num_rows = 3 in
459
+
let num_cols = 5 in
460
+
(* only first 5 columns are non-zero *)
461
+
let h_m = num_rows * m_sub in
462
+
(* 1536 *)
463
+
let h_n = num_cols * m_sub in
464
+
(* 2560 *)
465
+
(* Collect edges as (check_row, variable_col) pairs *)
466
+
let edges = Array.make h_m [] in
467
+
let add_edge chk var = edges.(chk) <- var :: edges.(chk) in
468
+
let process_layer proto =
469
+
for row_idx = 0 to num_rows - 1 do
470
+
for col_idx = 0 to num_cols - 1 do
471
+
let btype, bval = proto.(row_idx).(col_idx) in
472
+
if btype = hi then begin
473
+
(* Identity block shifted by bval *)
474
+
for c = 0 to m_sub - 1 do
475
+
let chk = (row_idx * m_sub) + c in
476
+
let var = (col_idx * m_sub) + ((c + bval) land mod_m) in
477
+
add_edge chk var
478
+
done
479
+
end
480
+
else if btype = hp then begin
481
+
(* Permutation block using theta/phi tables *)
482
+
for c = 0 to m_sub - 1 do
483
+
let j = c lsr log_circ in
484
+
(* j = c / circulant_size, 0..3 *)
485
+
let pi =
486
+
(((theta_k.(bval) + j) mod 4) lsl log_circ)
487
+
+ ((phi_j_k_m128.(j).(bval) + c) land mod_circ)
488
+
in
489
+
let chk = (row_idx * m_sub) + c in
490
+
let var = (col_idx * m_sub) + pi in
491
+
add_edge chk var
492
+
done
493
+
end
171
494
done
172
-
end
173
-
done;
495
+
done
496
+
in
497
+
process_layer tm_r12_proto_0;
498
+
process_layer tm_r12_proto_1;
499
+
process_layer tm_r12_proto_2;
500
+
(* Build sparse row representation *)
174
501
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))
502
+
Array.init h_m (fun i ->
503
+
let cols = List.sort_uniq compare edges.(i) in
504
+
Array.of_list cols)
178
505
in
179
-
let col_lists = Array.make n [] in
506
+
(* Build sparse column representation *)
507
+
let col_lists = Array.make h_n [] in
180
508
Array.iteri
181
509
(fun r cols ->
182
510
Array.iter (fun c -> col_lists.(c) <- r :: col_lists.(c)) cols)
···
184
512
let h_cols =
185
513
Array.map (fun rows -> Array.of_list (List.sort compare rows)) col_lists
186
514
in
187
-
(h_rows, h_cols)
515
+
(h_rows, h_cols, h_m, h_n)
188
516
189
-
let make_code ~n ~k ~m h_rows _h_cols =
190
-
let perm, gen_p = systematic_form ~m ~n h_rows in
517
+
(* Build a code from H matrix with column permutation for systematic form *)
518
+
let make_code_with_puncturing ~n_full ~k ~m ~punctured h_rows _h_cols =
519
+
let perm, gen_p = systematic_form ~m ~n:n_full h_rows in
191
520
(* Rebuild h_rows and h_cols with permuted columns *)
192
-
let inv_perm = Array.make n 0 in
521
+
let inv_perm = Array.make n_full 0 in
193
522
Array.iteri (fun new_pos orig -> inv_perm.(orig) <- new_pos) perm;
194
523
let h_rows' =
195
524
Array.map
···
199
528
cols')
200
529
h_rows
201
530
in
202
-
let col_lists = Array.make n [] in
531
+
let col_lists = Array.make n_full [] in
203
532
Array.iteri
204
533
(fun r cols ->
205
534
Array.iter (fun c -> col_lists.(c) <- r :: col_lists.(c)) cols)
···
207
536
let h_cols' =
208
537
Array.map (fun rows -> Array.of_list (List.sort compare rows)) col_lists
209
538
in
210
-
{ n; k; m; h_rows = h_rows'; h_cols = h_cols'; gen_p }
539
+
let n = n_full - punctured in
540
+
{ n; k; m; n_full; punctured; h_rows = h_rows'; h_cols = h_cols'; gen_p }
211
541
212
542
(* --- Rate 1/2 k=1024 code (lazy to avoid startup cost) --- *)
213
543
214
544
let ccsds_rate_1_2 =
215
545
let code =
216
546
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
547
+
let h_rows, h_cols, h_m, h_n = make_ar4ja_h () in
548
+
make_code_with_puncturing ~n_full:h_n ~k:1024 ~m:h_m ~punctured:512 h_rows
549
+
h_cols
224
550
end
225
551
in
226
552
Lazy.force code
···
257
583
Array.init code.m (fun j ->
258
584
Array.fold_left (fun acc i -> acc lxor info.(i)) 0 code.gen_p.(j))
259
585
in
260
-
(* Build codeword: [info | parity], packed into bytes *)
586
+
(* Build full codeword: [info | parity], n_full bits *)
587
+
let full_bits = Array.make code.n_full 0 in
588
+
Array.blit info 0 full_bits 0 code.k;
589
+
Array.blit parity 0 full_bits code.k code.m;
590
+
(* Puncture: remove the last 'punctured' bits to get n transmitted bits *)
261
591
let cw_bytes = (code.n + 7) / 8 in
262
592
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)
593
+
for i = 0 to code.n - 1 do
594
+
set_bit cw i full_bits.(i)
268
595
done;
269
596
cw
270
597
···
302
629
((code.n + 7) / 8)
303
630
cw_bits)
304
631
else begin
305
-
(* Initialize LLRs from hard-decision input *)
632
+
(* Initialize LLRs from hard-decision input.
633
+
For transmitted bits, use the channel LLR.
634
+
For punctured bits, use LLR = 0 (maximum uncertainty / erasure). *)
306
635
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)
636
+
Array.init code.n_full (fun i ->
637
+
if i < code.n then begin
638
+
(* Transmitted bit: read from codeword *)
639
+
if get_bit codeword i = 0 then hard_decision_llr
640
+
else Float.neg hard_decision_llr
641
+
end
642
+
else
643
+
(* Punctured bit: no channel information *)
644
+
0.0)
310
645
in
311
646
(* Message arrays.
312
647
For each edge (check i, variable j):
···
336
671
(* Build reverse index: for variable j and check i, what is the local
337
672
index of j in h_rows.(i)? *)
338
673
let var_to_check_local =
339
-
Array.init code.n (fun j ->
674
+
Array.init code.n_full (fun j ->
340
675
Array.map
341
676
(fun i ->
342
677
(* Find local index of j in h_rows.(i) *)
···
388
723
done;
389
724
390
725
(* --- Variable node update (vertical step) --- *)
391
-
for j = 0 to code.n - 1 do
726
+
for j = 0 to code.n_full - 1 do
392
727
let checks = var_to_check_local.(j) in
393
728
let num_checks_j = Array.length checks in
394
729
let total_r = ref channel_llr.(j) in
···
405
740
406
741
(* --- Hard decision and syndrome check --- *)
407
742
let hard =
408
-
Array.init code.n (fun j -> if posterior.(j) >= 0.0 then 0 else 1)
743
+
Array.init code.n_full (fun j -> if posterior.(j) >= 0.0 then 0 else 1)
409
744
in
410
745
if syndrome code hard then converged := true
411
746
done;
412
747
413
748
(* Extract the decoded information bits *)
414
749
let hard =
415
-
Array.init code.n (fun j -> if posterior.(j) >= 0.0 then 0 else 1)
750
+
Array.init code.n_full (fun j -> if posterior.(j) >= 0.0 then 0 else 1)
416
751
in
417
752
let data_bytes = (code.k + 7) / 8 in
418
753
let out = Bytes.make data_bytes '\x00' in
+16
-11
lib/ldpc.mli
+16
-11
lib/ldpc.mli
···
1
1
(** Low-Density Parity-Check codes with belief propagation decoding.
2
2
3
-
Generic LDPC encoder/decoder framework with systematic encoding and
4
-
sum-product (belief propagation) decoding. Configurable parity check matrix.
3
+
Implements the CCSDS 131.0-B-4 AR4JA (Accumulate-Repeat-4-Jagged-
4
+
Accumulate) codes with circulant-based parity-check matrices, systematic
5
+
encoding, and sum-product belief propagation decoding.
5
6
6
7
@see <https://public.ccsds.org/Pubs/131x0b4.pdf> CCSDS 131.0-B-4 Section 7
7
8
@see <https://public.ccsds.org/Pubs/231x0b4.pdf> CCSDS 231.0-B Section 4 *)
···
11
12
parameters (n, k, m). *)
12
13
13
14
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.
15
+
(** Rate 1/2 LDPC code with k=1024, n=2048.
16
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
-
*)
17
+
Implements the CCSDS 131.0-B-4 AR4JA (Accumulate-Repeat-4-Jagged-
18
+
Accumulate) code. The parity-check matrix is constructed from circulant
19
+
permutation blocks as specified in the standard, using a two-step expansion
20
+
of the AR4JA protograph.
21
+
22
+
The internal code has n_full=2560 variable nodes and m=1536 check equations.
23
+
512 bits are punctured (not transmitted), giving the effective rate 1/2 code
24
+
with k=1024 information bits and n=2048 transmitted bits. *)
22
25
23
26
val encode : code -> bytes -> bytes
24
27
(** [encode code data] performs systematic LDPC encoding.
25
28
26
29
The input [data] must contain at least [k] bits (i.e., [ceil(k/8)] bytes).
27
30
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.
31
+
bits followed by parity bits, packed MSB-first. For codes with punctured
32
+
bits, the punctured positions are removed from the output.
29
33
30
34
@raise Invalid_argument if the input is too short. *)
31
35
···
34
38
sum-product (belief propagation) algorithm.
35
39
36
40
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
41
+
converted to LLRs (bit 0 -> +LLR, bit 1 -> -LLR). Punctured bit positions
42
+
are initialized with LLR=0 (maximum uncertainty). The decoder iterates until
38
43
the syndrome is zero or [max_iter] iterations are reached (default 50).
39
44
40
45
Returns [Ok data] with the [ceil(k/8)] decoded information bytes on success,