comnum_fec/ccsds/
ldpc.rs

1//! The LDPC codes of CCSDS 131.0-B-5 Sec. 7: the (8160, 7136) code of near-Earth missions
2//! and the AR4JA codes of deep-space missions, [`CcsdsLdpc`].
3
4use crate::block::BlockCode;
5use crate::block::SoftDecode;
6use crate::block::container_traits;
7use crate::block::expect_len;
8use crate::ldpc::LdpcCode;
9use comnum_core::Error;
10use comnum_core::Result;
11use comnum_core::ensure;
12use comnum_core::signal::KNOWN_BIT_LLR;
13use comnum_core::signal::hard_decisions;
14use comnum_math::gf2;
15
16/// Table 7-1: the two positions of the ones in the first row of every circulant $A_{i,j}$
17/// of the (8176, 7156) code, $i = 1, 2$, $j = 1, \dots, 16$.
18const LDPC_C2_A: [[[u16; 2]; 16]; 2] = [
19    [
20        [0, 176],
21        [12, 239],
22        [0, 352],
23        [24, 431],
24        [0, 392],
25        [151, 409],
26        [0, 351],
27        [9, 359],
28        [0, 307],
29        [53, 329],
30        [0, 207],
31        [18, 281],
32        [0, 399],
33        [202, 457],
34        [0, 247],
35        [36, 261],
36    ],
37    [
38        [99, 471],
39        [130, 473],
40        [198, 435],
41        [260, 478],
42        [215, 420],
43        [282, 481],
44        [48, 396],
45        [193, 445],
46        [273, 430],
47        [302, 451],
48        [96, 379],
49        [191, 386],
50        [244, 467],
51        [364, 470],
52        [51, 382],
53        [192, 414],
54    ],
55];
56
57/// Table C-1: the first rows of the circulants $B_{i,j}$ of the generator of the
58/// (8176, 7154) subcode, $i = 1, \dots, 14$, $j = 1, 2$: 128 hexadecimal digits each, the
59/// first one octal (3 bits), 511 bits, position 0 first.
60const LDPC_C2_B: [[&str; 2]; 14] = [
61    [
62        "55BF56CC55283DFEEFEA8C8CFF04E1EBD9067710988E25048D67525426939E2068D2DC6FCD2F822BEB6BD96C8A76F4932AAE9BC53AD20A2A9C86BB461E43759C",
63        "6855AE08698A50AA3051768793DC238544AF3FE987391021AAF6383A6503409C3CE971A80B3ECE12363EE809A01D91204F1811123EAB867D3E40E8C652585D28",
64    ],
65    [
66        "62B21CF0AEE0649FA67B7D0EA6551C1CD194CA77501E0FCF8C85867B9CF679C18BCF7939E10F8550661848A4E0A9E9EDB7DAB9EDABA18C168C8E28AACDDEAB1E",
67        "64B71F486AD57125660C4512247B229F0017BA649C6C11148FB00B70808286F1A9790748D296A593FA4FD2C6D7AAF7750F0C71B31AEE5B400C7F5D73AAF00710",
68    ],
69    [
70        "681A8E51420BD8294ECE13E491D618083FFBBA830DB5FAF330209877D801F92B5E07117C57E75F6F0D873B3E520F21EAFD78C1612C6228111A369D5790F5929A",
71        "04DF1DD77F1C20C1FB570D7DD7A1219EAECEA4B2877282651B0FFE713DF338A63263BC0E324A87E2DC1AD64C9F10AAA585ED6905946EE167A73CF04AD2AF9218",
72    ],
73    [
74        "35951FEE6F20C902296C9488003345E6C5526C5519230454C556B8A04FC0DC642D682D94B4594B5197037DF15B5817B26F16D0A3302C09383412822F6D2B234E",
75        "7681CF7F278380E28F1262B22F40BF3405BFB92311A8A34D084C086464777431DBFDDD2E82A2E6742BAD6533B51B2BDEE0377E9F6E63DCA0B0F1DF97E73D5CD8",
76    ],
77    [
78        "188157AE41830744BAE0ADA6295E08B79A44081E111F69BBE7831D07BEEBF76232E065F752D4F218D39B6C5BF20AE5B8FF172A7F1F680E6BF5AAC3C4343736C2",
79        "5D80A6007C175B5C0DD88A442440E2C29C6A136BBCE0D95A58A83B48CA0E7474E9476C92E33D164BFF943A61CE1031DFF441B0B175209B498394F4794644392E",
80    ],
81    [
82        "60CD1F1C282A1612657E8C7C1420332CA245C0756F78744C807966C3E1326438878BD2CCC83388415A612705AB192B3512EEF0D95248F7B73E5B0F412BF76DB4",
83        "434B697B98C9F3E48502C8DBD891D0A0386996146DEBEF11D4B833033E05EDC28F808F25E8F314135E6675B7608B66F7FF3392308242930025DDC4BB65CD7B6E",
84    ],
85    [
86        "766855125CFDC804DAF8DBE3660E8686420230ED4E049DF11D82E357C54FE256EA01F5681D95544C7A1E32B7C30A8E6CF5D0869E754FFDE6AEFA6D7BE8F1B148",
87        "222975D325A487FE560A6D146311578D9C5501D28BC0A1FB48C9BDA173E869133A3AA9506C42AE9F466E85611FC5F8F74E439638D66D2F00C682987A96D8887C",
88    ],
89    [
90        "14B5F98E8D55FC8E9B4EE453C6963E052147A857AC1E08675D99A308E7269FAC5600D7B155DE8CB1BAC786F45B46B523073692DE745FDF10724DDA38FD093B1C",
91        "1B71AFFB8117BCF8B5D002A99FEEA49503C0359B056963FE5271140E626F6F8FCE9F29B37047F9CA89EBCE760405C6277F329065DF21AB3B779AB3E8C8955400",
92    ],
93    [
94        "0008B4E899E5F7E692BDCE69CE3FAD997183CFAEB2785D0C3D9CAE510316D4BD65A2A06CBA7F4E4C4A80839ACA81012343648EEA8DBBA2464A68E115AB3F4034",
95        "5B7FE6808A10EA42FEF0ED9B41920F82023085C106FBBC1F56B567A14257021BC5FDA60CBA05B08FAD6DC3B0410295884C7CCDE0E56347D649DE6DDCEEB0C95E",
96    ],
97    [
98        "5E9B2B33EF82D0E64AA2226D6A0ADCD179D5932EE1CF401B336449D0FF775754CA56650716E61A43F963D59865C7F017F53830514306649822CAA72C152F6EB2",
99        "2CD8140C8A37DE0D0261259F63AA2A420A8F81FECB661DBA5C62DF6C817B4A61D2BC1F068A50DFD0EA8FE1BD387601062E2276A4987A19A70B460C54F215E184",
100    ],
101    [
102        "06F1FF249192F2EAF063488E267EEE994E7760995C4FA6FFA0E4241825A7F5B65C74FB16AC4C891BC008D33AD4FF97523EE5BD14126916E0502FF2F8E4A07FC2",
103        "65287840D00243278F41CE1156D1868F24E02F91D3A1886ACE906CE741662B40B4EFDFB90F76C1ADD884D920AFA8B3427EEB84A759FA02E00635743F50B942F0",
104    ],
105    [
106        "4109DA2A24E41B1F375645229981D4B7E88C36A12DAB64E91C764CC43CCEC188EC8C5855C8FF488BB91003602BEF43DBEC4A621048906A2CDC5DBD4103431DB8",
107        "2185E3BC7076BA51AAD6B199C8C60BCD70E8245B874927136E6D8DD527DF0693DC10A1C8E51B5BE93FF7538FA138B335738F4315361ABF8C73BF40593AE22BE4",
108    ],
109    [
110        "228845775A262505B47288E065B23B4A6D78AFBDDB2356B392C692EF56A35AB4AA27767DE72F058C6484457C95A8CCDD0EF225ABA56B7657B7F0E947DC17F972",
111        "2630C6F79878E50CF5ABD353A6ED80BEACC7169179EA57435E44411BC7D566136DFA983019F3443DE8E4C60940BC4E31DCEAD514D755AF95A622585D69572692",
112    ],
113    [
114        "7273E8342918E097B1C1F5FEF32A150AEF5E11184782B5BD5A1D8071E94578B0AC722D7BF49E8C78D391294371FFBA7B88FABF8CC03A62B940CE60D669DFB7B6",
115        "087EA12042793307045B283D7305E93D8F74725034E77D25D3FF043ADC5F8B5B186DB70A968A816835EFB575952EAE7EA4E76DF0D5F097590E1A2A978025573E",
116    ],
117];
118
119/// $\theta_{k}$ of Tables 7-3 and 7-4, $k = 1, \dots, 26$.
120const AR4JA_THETA: [u8; 26] = [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];
121
122/// $\phi_{k}(j, M)$ of Tables 7-3 ($j = 0, 1$) and 7-4 ($j = 2, 3$), for $M = 2^{7}, \dots,
123/// 2^{13}$.
124const AR4JA_PHI: [[[u16; 7]; 4]; 26] = [
125    [[1, 59, 16, 160, 108, 226, 1148], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0]],
126    [
127        [22, 18, 103, 241, 126, 618, 2032],
128        [27, 32, 53, 182, 375, 767, 1822],
129        [12, 46, 8, 35, 219, 254, 318],
130        [13, 44, 35, 162, 312, 285, 1189],
131    ],
132    [
133        [0, 52, 105, 185, 238, 404, 249],
134        [30, 21, 74, 249, 436, 227, 203],
135        [30, 45, 119, 167, 16, 790, 494],
136        [19, 51, 97, 7, 503, 554, 458],
137    ],
138    [
139        [26, 23, 0, 251, 481, 32, 1807],
140        [28, 36, 45, 65, 350, 247, 882],
141        [18, 27, 89, 214, 263, 642, 1467],
142        [14, 12, 112, 31, 388, 809, 460],
143    ],
144    [
145        [0, 11, 50, 209, 96, 912, 485],
146        [7, 30, 47, 70, 260, 284, 1989],
147        [10, 48, 31, 84, 415, 248, 757],
148        [15, 15, 64, 164, 48, 185, 1039],
149    ],
150    [
151        [10, 7, 29, 103, 28, 950, 1044],
152        [1, 29, 0, 141, 84, 370, 957],
153        [16, 37, 122, 206, 403, 899, 1085],
154        [20, 12, 93, 11, 7, 49, 1000],
155    ],
156    [
157        [5, 22, 115, 90, 59, 534, 717],
158        [8, 44, 59, 237, 318, 482, 1705],
159        [13, 41, 1, 122, 184, 328, 1630],
160        [17, 4, 99, 237, 185, 101, 1265],
161    ],
162    [
163        [18, 25, 30, 184, 225, 63, 873],
164        [20, 29, 102, 77, 382, 273, 1083],
165        [9, 13, 69, 67, 279, 518, 64],
166        [4, 7, 94, 125, 328, 82, 1223],
167    ],
168    [
169        [3, 27, 92, 248, 323, 971, 364],
170        [26, 39, 25, 55, 169, 886, 1072],
171        [7, 9, 92, 147, 198, 477, 689],
172        [4, 2, 103, 133, 254, 898, 874],
173    ],
174    [
175        [22, 30, 78, 12, 28, 304, 1926],
176        [24, 14, 3, 12, 213, 634, 354],
177        [15, 49, 47, 54, 307, 404, 1300],
178        [11, 30, 91, 99, 202, 627, 1292],
179    ],
180    [
181        [3, 43, 70, 111, 386, 409, 1241],
182        [4, 22, 88, 227, 67, 762, 1942],
183        [16, 36, 11, 23, 432, 698, 148],
184        [17, 53, 3, 105, 285, 154, 1491],
185    ],
186    [
187        [8, 14, 66, 66, 305, 708, 1769],
188        [12, 15, 65, 42, 313, 184, 446],
189        [18, 10, 31, 93, 240, 160, 777],
190        [20, 23, 6, 17, 11, 65, 631],
191    ],
192    [
193        [25, 46, 39, 173, 34, 719, 532],
194        [23, 48, 62, 52, 242, 696, 1456],
195        [4, 11, 19, 20, 454, 497, 1431],
196        [8, 29, 39, 97, 168, 81, 464],
197    ],
198    [
199        [25, 62, 84, 42, 510, 176, 768],
200        [15, 55, 68, 243, 188, 413, 1940],
201        [23, 18, 66, 197, 294, 100, 659],
202        [22, 37, 113, 91, 127, 823, 461],
203    ],
204    [
205        [2, 44, 79, 157, 147, 743, 1138],
206        [15, 39, 91, 179, 1, 854, 1660],
207        [5, 54, 49, 46, 479, 518, 352],
208        [19, 42, 92, 211, 8, 50, 844],
209    ],
210    [
211        [27, 12, 70, 174, 199, 759, 965],
212        [22, 11, 70, 250, 306, 544, 1661],
213        [3, 40, 81, 162, 289, 92, 1177],
214        [15, 48, 119, 128, 437, 413, 392],
215    ],
216    [
217        [7, 38, 29, 104, 347, 674, 141],
218        [31, 1, 115, 247, 397, 864, 587],
219        [29, 27, 96, 101, 373, 464, 836],
220        [5, 4, 74, 82, 475, 462, 922],
221    ],
222    [
223        [7, 47, 32, 144, 391, 958, 1527],
224        [3, 50, 31, 164, 80, 82, 708],
225        [11, 35, 38, 76, 104, 592, 1572],
226        [21, 10, 73, 115, 85, 175, 256],
227    ],
228    [
229        [15, 1, 45, 43, 165, 984, 505],
230        [29, 40, 121, 17, 33, 1009, 1466],
231        [4, 25, 83, 78, 141, 198, 348],
232        [17, 18, 116, 248, 419, 715, 1986],
233    ],
234    [
235        [10, 52, 113, 181, 414, 11, 1312],
236        [21, 62, 45, 31, 7, 437, 433],
237        [8, 46, 42, 253, 270, 856, 1040],
238        [9, 56, 31, 62, 459, 537, 19],
239    ],
240    [
241        [4, 61, 86, 250, 97, 413, 1840],
242        [2, 27, 56, 149, 447, 36, 1345],
243        [2, 24, 58, 124, 439, 235, 779],
244        [20, 9, 127, 26, 468, 722, 266],
245    ],
246    [
247        [19, 10, 1, 202, 158, 925, 709],
248        [5, 38, 54, 105, 336, 562, 867],
249        [11, 33, 24, 143, 333, 134, 476],
250        [18, 11, 98, 140, 209, 37, 471],
251    ],
252    [
253        [7, 55, 42, 68, 86, 687, 1427],
254        [11, 40, 108, 183, 424, 816, 1551],
255        [11, 18, 25, 63, 399, 542, 191],
256        [31, 23, 23, 121, 311, 488, 1166],
257    ],
258    [
259        [9, 7, 118, 177, 168, 752, 989],
260        [26, 15, 14, 153, 134, 452, 2041],
261        [3, 37, 92, 41, 14, 545, 1393],
262        [13, 8, 38, 12, 211, 179, 1300],
263    ],
264    [
265        [26, 12, 33, 170, 506, 867, 1925],
266        [9, 11, 30, 177, 152, 290, 1383],
267        [15, 35, 38, 214, 277, 777, 1752],
268        [2, 7, 18, 41, 510, 430, 1033],
269    ],
270    [
271        [17, 2, 126, 89, 489, 323, 270],
272        [17, 18, 116, 19, 492, 778, 1790],
273        [13, 21, 120, 70, 412, 483, 1627],
274        [18, 24, 62, 249, 320, 264, 1606],
275    ],
276];
277
278/// The size of a circulant of the (8176, 7156) code.
279const CIRCULANT: usize = 511;
280
281/// The code rates of the AR4JA codes (Sec. 7.4).
282#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
283pub enum Ar4jaRate {
284    /// 1/2, $K = 2$.
285    R1_2,
286    /// 2/3, $K = 4$.
287    R2_3,
288    /// 4/5, $K = 8$.
289    R4_5,
290}
291
292impl Ar4jaRate {
293    /// The number $K$ of information block columns of $H$.
294    pub fn blocks(self) -> usize {
295        match self {
296            Ar4jaRate::R1_2 => 2,
297            Ar4jaRate::R2_3 => 4,
298            Ar4jaRate::R4_5 => 8,
299        }
300    }
301
302    /// The name of the rate, `1/2`, `2/3` or `4/5`.
303    pub fn name(self) -> &'static str {
304        match self {
305            Ar4jaRate::R1_2 => "1/2",
306            Ar4jaRate::R2_3 => "2/3",
307            Ar4jaRate::R4_5 => "4/5",
308        }
309    }
310}
311
312impl std::str::FromStr for Ar4jaRate {
313    type Err = Error;
314
315    /// `1/2`, `2/3` or `4/5`.
316    fn from_str(s: &str) -> Result<Ar4jaRate> {
317        match s.trim() {
318            "1/2" => Ok(Ar4jaRate::R1_2),
319            "2/3" => Ok(Ar4jaRate::R2_3),
320            "4/5" => Ok(Ar4jaRate::R4_5),
321            _ => Err(Error::Config(format!("the AR4JA rate must be 1/2, 2/3 or 4/5, got {s:?}"))),
322        }
323    }
324}
325
326/// The permutation $\pi_{k}(i)$ of $\Pi_{k}$, $k = 1, \dots, 26$, for the submatrix size $M$
327/// ($2^{7}$ to $2^{13}$), Sec. 7.4.2.4.
328pub(super) fn ar4ja_pi(k: usize, i: usize, m: usize) -> usize {
329    let col = m.trailing_zeros() as usize - 7;
330    let q = m / 4;
331    let j = 4 * i / m;
332    q * ((usize::from(AR4JA_THETA[k - 1]) + j) % 4) + (usize::from(AR4JA_PHI[k - 1][j][col]) + i) % q
333}
334
335/// A block of $H$ of an AR4JA code: the sum of an identity (if `identity`) and of the
336/// permutations $\Pi_{k}$ of `pis`.
337#[derive(Clone, Copy)]
338struct Block {
339    identity: bool,
340    pis: &'static [usize],
341}
342
343const ZERO: Block = Block { identity: false, pis: &[] };
344
345const IDENTITY: Block = Block { identity: true, pis: &[] };
346
347const fn pis(pis: &'static [usize]) -> Block {
348    Block { identity: false, pis }
349}
350
351/// The block rows of $H$ of an AR4JA code, $3 \times (K + 3)$ blocks (Secs. 7.4.2.2,
352/// 7.4.2.3).
353fn ar4ja_blocks(rate: Ar4jaRate) -> Vec<[Block; 3]> {
354    let h12 = [
355        [ZERO, IDENTITY, IDENTITY],
356        [ZERO, IDENTITY, pis(&[5, 6])],
357        [IDENTITY, ZERO, ZERO],
358        [ZERO, IDENTITY, pis(&[7, 8])],
359        [Block { identity: true, pis: &[1] }, pis(&[2, 3, 4]), IDENTITY],
360    ];
361    let p23 = [[ZERO, pis(&[9, 10, 11]), IDENTITY], [ZERO, IDENTITY, pis(&[12, 13, 14])]];
362    let p45 = [
363        [ZERO, pis(&[21, 22, 23]), IDENTITY],
364        [ZERO, IDENTITY, pis(&[24, 25, 26])],
365        [ZERO, pis(&[15, 16, 17]), IDENTITY],
366        [ZERO, IDENTITY, pis(&[18, 19, 20])],
367    ];
368    // block columns, each the three blocks of its block rows
369    let mut cols: Vec<[Block; 3]> = Vec::new();
370    if rate == Ar4jaRate::R4_5 {
371        cols.extend(p45);
372    }
373    if rate != Ar4jaRate::R1_2 {
374        cols.extend(p23);
375    }
376    cols.extend(h12);
377    cols
378}
379
380/// The LDPC codes of CCSDS (see the module): the decoder of [`LdpcCode`] on the full
381/// parity-check matrix, the encoder of the standard, the virtual zeros and the punctured
382/// bits.
383#[derive(Clone, Debug, PartialEq)]
384pub struct CcsdsLdpc {
385    /// The decoder: the code of the full parity-check matrix.
386    code: LdpcCode,
387    kind: LdpcKind,
388}
389
390#[derive(Clone, Debug, PartialEq)]
391enum LdpcKind {
392    /// The (8160, 7136) code: the 28 generator circulants of Table C-1 (511 bits each).
393    C2 { b: Vec<Vec<u8>> },
394    /// An AR4JA code: $M$, $K$, the permutations of the blocks and the inverse of
395    /// $C = I + (\Pi_{7} + \Pi_{8})(\Pi_{2} + \Pi_{3} + \Pi_{4})$.
396    Ar4ja { m: usize, blocks: usize, perms: Vec<Vec<usize>>, c_inv: gf2::Matrix },
397}
398
399impl CcsdsLdpc {
400    /// The (8160, 7136) code of rate 223/255 (Sec. 7.3).
401    ///
402    /// # Errors
403    /// Never: the tables are those of the standard.
404    pub fn c2() -> Result<CcsdsLdpc> {
405        let mut rows = Vec::with_capacity(2 * CIRCULANT);
406        for a in &LDPC_C2_A {
407            for r in 0..CIRCULANT {
408                rows.push(
409                    a.iter()
410                        .enumerate()
411                        .flat_map(|(j, pos)| pos.iter().map(move |&p| CIRCULANT * j + (usize::from(p) + r) % CIRCULANT))
412                        .collect(),
413                );
414            }
415        }
416        let code = LdpcCode::without_encoder(16 * CIRCULANT, rows, (0..14 * CIRCULANT).collect())?;
417        let b = LDPC_C2_B.iter().flatten().map(|h| circulant_row(h)).collect::<Result<_>>()?;
418        Ok(CcsdsLdpc { code, kind: LdpcKind::C2 { b } })
419    }
420
421    /// The AR4JA code of `k` information bits (1024, 4096 or 16384) and rate `rate`
422    /// (Sec. 7.4). Building the code of $k = 16384$ and rate 1/2 inverts an $8192 \times
423    /// 8192$ binary matrix (seconds in release).
424    ///
425    /// # Errors
426    /// [`Error::Config`] for another $k$.
427    pub fn ar4ja(k: usize, rate: Ar4jaRate) -> Result<CcsdsLdpc> {
428        ensure(matches!(k, 1024 | 4096 | 16384), || format!("the AR4JA codes have k = 1024, 4096 or 16384, got {k}"))?;
429        let kb = rate.blocks();
430        let m = k / kb;
431        let perms: Vec<Vec<usize>> = (1..=26).map(|p| (0..m).map(|i| ar4ja_pi(p, i, m)).collect()).collect();
432        let cols = ar4ja_blocks(rate);
433        let n_full = cols.len() * m;
434        let mut rows: Vec<Vec<usize>> = vec![Vec::new(); 3 * m];
435        for (c, col) in cols.iter().enumerate() {
436            for (r, block) in col.iter().enumerate() {
437                for i in 0..m {
438                    let row = &mut rows[r * m + i];
439                    if block.identity {
440                        toggle(row, c * m + i);
441                    }
442                    for &p in block.pis {
443                        toggle(row, c * m + perms[p - 1][i]);
444                    }
445                }
446            }
447        }
448        for row in &mut rows {
449            row.sort_unstable();
450        }
451        // C = I + (P7 + P8)(P2 + P3 + P4): its row i lists i and pi_a(pi_b(i))
452        let c_rows: Vec<Vec<usize>> = (0..m)
453            .map(|i| {
454                let mut row = vec![i];
455                for &b in &[7, 8] {
456                    let j = perms[b - 1][i];
457                    row.extend([2, 3, 4].iter().map(|&a| perms[a - 1][j]));
458                }
459                row
460            })
461            .collect();
462        let c_inv = gf2::Matrix::from_rows(m, &c_rows)
463            .inverse()
464            .ok_or_else(|| Error::Numerical("the AR4JA matrix C is singular".into()))?;
465        let code = LdpcCode::without_encoder(n_full, rows, (0..kb * m).collect())?;
466        Ok(CcsdsLdpc { code, kind: LdpcKind::Ar4ja { m, blocks: kb, perms, c_inv } })
467    }
468
469    /// The code of the full parity-check matrix (its decoder settings are public fields).
470    pub fn ldpc(&self) -> &LdpcCode {
471        &self.code
472    }
473
474    /// The code of the full parity-check matrix, to change its decoder settings.
475    pub fn ldpc_mut(&mut self) -> &mut LdpcCode {
476        &mut self.code
477    }
478
479    /// The codeword of the full parity-check matrix ($n$ = 8176 or $(K + 3) M$ bits) of
480    /// `info`: the 18 virtual zeros of the rate-7/8 code, the information, the parity
481    /// bits, the punctured ones of the AR4JA codes included; [`LdpcCode::is_codeword`] of
482    /// [`CcsdsLdpc::ldpc`] checks it.
483    ///
484    /// # Errors
485    /// [`Error::Shape`] unless `info` has $k$ bits.
486    pub fn full_codeword(&self, info: &[u8]) -> Result<Vec<u8>> {
487        expect_len(info.len(), self.k(), "bits")?;
488        Ok(self.full_codeword_unchecked(info))
489    }
490
491    fn full_codeword_unchecked(&self, info: &[u8]) -> Vec<u8> {
492        match &self.kind {
493            LdpcKind::C2 { b } => {
494                let mut c = vec![0u8; 18];
495                c.extend_from_slice(info);
496                let mut parity = vec![0u8; 2 * CIRCULANT];
497                for i in 0..14 {
498                    for r in 0..CIRCULANT {
499                        if c[i * CIRCULANT + r] == 1 {
500                            for j in 0..2 {
501                                let row = &b[2 * i + j];
502                                for (t, p) in parity[j * CIRCULANT..(j + 1) * CIRCULANT].iter_mut().enumerate() {
503                                    *p ^= row[(t + CIRCULANT - r) % CIRCULANT];
504                                }
505                            }
506                        }
507                    }
508                }
509                c.extend(parity);
510                c
511            }
512            LdpcKind::Ar4ja { m, blocks, perms, c_inv } => {
513                let (m, kb) = (*m, *blocks);
514                // s = Q u: the syndrome of the information block columns
515                let mut s = vec![0u8; 3 * m];
516                for (r, row) in self.code.rows().iter().enumerate() {
517                    s[r] = row.iter().filter(|&&v| v < kb * m).fold(0, |a, &v| a ^ info[v]);
518                }
519                let pi = |p: usize, v: &[u8]| -> Vec<u8> { (0..m).map(|i| v[perms[p - 1][i]]).collect() };
520                let xor = |a: &mut [u8], b: &[u8]| a.iter_mut().zip(b).for_each(|(x, y)| *x ^= y);
521                let (s0, s1, s2) = (&s[..m], &s[m..2 * m], &s[2 * m..]);
522                // x4 = C^-1 (s2 + (P7 + P8) s1)
523                let mut y = s2.to_vec();
524                xor(&mut y, &pi(7, s1));
525                xor(&mut y, &pi(8, s1));
526                let x4 = c_inv.mul_vec(&y);
527                // x3 = s1 + (P2 + P3 + P4) x4, x2 = s0 + (I + P1) x4
528                let mut x3 = s1.to_vec();
529                for p in [2, 3, 4] {
530                    xor(&mut x3, &pi(p, &x4));
531                }
532                let mut x2 = s0.to_vec();
533                xor(&mut x2, &x4);
534                xor(&mut x2, &pi(1, &x4));
535                let mut c = info.to_vec();
536                c.extend(x2);
537                c.extend(x3);
538                c.extend(x4);
539                c
540            }
541        }
542    }
543
544    /// The channel LLRs of the full codeword from the received ones: the virtual zeros
545    /// known, the punctured bits at zero.
546    fn full_llr(&self, llr: &[f64]) -> Vec<f64> {
547        match &self.kind {
548            LdpcKind::C2 { .. } => {
549                let mut full = vec![KNOWN_BIT_LLR; 18];
550                full.extend_from_slice(&llr[..16 * CIRCULANT - 18]);
551                full
552            }
553            LdpcKind::Ar4ja { m, .. } => {
554                let mut full = llr.to_vec();
555                full.resize(llr.len() + m, 0.0);
556                full
557            }
558        }
559    }
560
561    /// The information positions in the full codeword.
562    fn info_range(&self) -> std::ops::Range<usize> {
563        match &self.kind {
564            LdpcKind::C2 { .. } => 18..18 + 7136,
565            LdpcKind::Ar4ja { m, blocks, .. } => 0..m * blocks,
566        }
567    }
568}
569
570/// Toggle `v` in a row of $H$ (two equal terms of a block cancel).
571fn toggle(row: &mut Vec<usize>, v: usize) {
572    match row.iter().position(|&x| x == v) {
573        Some(p) => {
574            row.swap_remove(p);
575        }
576        None => row.push(v),
577    }
578}
579
580/// The 511 bits of a first row of Table C-1: the 128 hexadecimal digits, the first one
581/// read as 3 bits.
582fn circulant_row(hex: &str) -> Result<Vec<u8>> {
583    let mut bits = Vec::with_capacity(512);
584    for (i, ch) in hex.chars().enumerate() {
585        let d = ch.to_digit(16).ok_or_else(|| Error::Config(format!("not a hexadecimal digit: {ch:?}")))?;
586        let width = if i == 0 { 3 } else { 4 };
587        bits.extend((0..width).rev().map(|k| ((d >> k) & 1) as u8));
588    }
589    ensure(bits.len() == CIRCULANT, || format!("a circulant row of {} bits", bits.len()))?;
590    Ok(bits)
591}
592
593impl BlockCode for CcsdsLdpc {
594    fn k(&self) -> usize {
595        self.info_range().len()
596    }
597
598    /// The transmitted bits: 8160, or $(K + 2) M$ (Table 7-5).
599    fn n(&self) -> usize {
600        match &self.kind {
601            LdpcKind::C2 { .. } => 8160,
602            LdpcKind::Ar4ja { m, blocks, .. } => (blocks + 2) * m,
603        }
604    }
605
606    fn encode_block(&self, info: &[u8]) -> Result<Vec<u8>> {
607        expect_len(info.len(), self.k(), "bits")?;
608        let mut c = self.full_codeword_unchecked(info);
609        match &self.kind {
610            LdpcKind::C2 { .. } => {
611                c.drain(..18);
612                c.extend([0, 0]);
613            }
614            LdpcKind::Ar4ja { m, .. } => c.truncate(c.len() - m),
615        }
616        Ok(c)
617    }
618
619    fn decode_block(&self, llr: &[f64]) -> Result<Vec<u8>> {
620        Ok(hard_decisions(&self.app(llr, None)?))
621    }
622}
623
624impl SoftDecode for CcsdsLdpc {
625    /// The a posteriori LLRs of the information bits by belief propagation on the full
626    /// parity-check matrix; the a priori LLRs add to the channel LLRs of the information
627    /// bits.
628    fn app_into(&self, llr: &[f64], apriori: Option<&[f64]>, out: &mut [f64]) -> Result<usize> {
629        expect_len(llr.len(), self.n(), "LLRs")?;
630        expect_len(out.len(), self.k(), "a posteriori LLRs")?;
631        let mut full = self.full_llr(llr);
632        let range = self.info_range();
633        if let Some(a) = apriori {
634            expect_len(a.len(), self.k(), "a priori LLRs")?;
635            for (l, a) in full[range.clone()].iter_mut().zip(a) {
636                *l += a;
637            }
638        }
639        let (app, iterations) = self.code.decode_llr(&full)?;
640        out.copy_from_slice(&app[range]);
641        Ok(iterations)
642    }
643}
644
645container_traits!(CcsdsLdpc);
comnum_fec/ccsds/ldpc.rs

comnum_fec/ccsds/
ldpc.rs
