ldpc_toolbox/simulation/
ber.rs

1//! BER simulation
2//!
3//! This module contains utilities for BER simulation.
4
5use super::{
6    channel::{AwgnChannel, Channel},
7    factory::Ber,
8    interleaving::Interleaver,
9    modulation::{Demodulator, Modulation, Modulator},
10    puncturing::Puncturer,
11};
12use crate::{
13    decoder::{
14        LdpcDecoder,
15        factory::{DecoderFactory, DecoderImplementation},
16    },
17    encoder::{Encoder, Error},
18    gf2::GF2,
19    sparse::SparseMatrix,
20};
21use ndarray::Array1;
22use num_traits::{One, Zero};
23use rand::{Rng, distr::StandardUniform};
24use std::{
25    sync::mpsc::{self, Receiver, Sender, SyncSender, TryRecvError},
26    time::{Duration, Instant},
27};
28
29/// BER test.
30///
31/// This struct is used to configure and run a BER test.
32#[derive(Debug)]
33pub struct BerTest<Mod: Modulation, Dec = DecoderImplementation> {
34    decoder_implementation: Dec,
35    h: SparseMatrix,
36    num_workers: usize,
37    k: usize,
38    n: usize,
39    n_cw: usize,
40    rate: f64,
41    encoder: Encoder,
42    puncturer: Option<Puncturer>,
43    interleaver: Option<Interleaver>,
44    modulator: Mod::Modulator,
45    ebn0s_db: Vec<f32>,
46    statistics: Vec<Statistics>,
47    bch_max_errors: u64,
48    max_iterations: usize,
49    max_frame_errors: u64,
50    min_run_time: Duration,
51    max_run_time: Duration,
52    reporter: Option<Reporter>,
53    last_reported: Instant,
54}
55
56/// BER test parameters.
57///
58/// This struct contains the configuration of a BER test.
59#[derive(Debug)]
60pub struct BerTestParameters<'a, Dec = DecoderImplementation> {
61    /// LDPC parity check matrix.
62    pub h: SparseMatrix,
63    /// LDPC decoder implementation.
64    pub decoder_implementation: Dec,
65    /// Codeword puncturing pattern.
66    pub puncturing_pattern: Option<&'a [bool]>,
67    /// Codeword interleaving.
68    ///
69    /// A negative value indicates that the columns should be read backwards.
70    pub interleaving_columns: Option<isize>,
71    /// Maximum number of frame errors per Eb/N0.
72    pub max_frame_errors: u64,
73    /// Minimum run time per Eb/N0.
74    ///
75    /// If this value is set to `None`, the simulation finishes as soon as the
76    /// frame errors are reached.
77    pub min_run_time: Option<Duration>,
78    /// Maximum run time per Eb/N0.
79    ///
80    /// If this value is set to `None`, the simulation only finishes when the
81    /// frame errors are reached.
82    pub max_run_time: Option<Duration>,
83    /// Maximum number of iterations per codeword.
84    pub max_iterations: usize,
85    /// List of Eb/N0's (in dB) to simulate.
86    pub ebn0s_db: &'a [f32],
87    /// An optional reporter object to which the BER test will send periodic
88    /// updates about its progress.
89    pub reporter: Option<Reporter>,
90    /// Maximum number of bit errors that the BCH decoder can correct.
91    ///
92    /// A value of zero means that there is no BCH decoder.
93    pub bch_max_errors: u64,
94    /// Number of worker threads.
95    pub num_workers: usize,
96}
97
98#[derive(Debug)]
99struct Worker<Mod: Modulation> {
100    terminate_rx: Receiver<()>,
101    results_tx: Sender<WorkerResult>,
102    k: usize,
103    encoder: Encoder,
104    puncturer: Option<Puncturer>,
105    interleaver: Option<Interleaver>,
106    modulator: Mod::Modulator,
107    channel: AwgnChannel,
108    demodulator: Mod::Demodulator,
109    decoder: Box<dyn LdpcDecoder>,
110    max_iterations: usize,
111}
112
113#[derive(Debug, Clone)]
114struct WorkerResultOk {
115    bit_errors: u64,
116    frame_error: bool,
117    false_decode: bool,
118    iterations: u64,
119}
120
121type WorkerResult = Result<WorkerResultOk, ()>;
122
123#[derive(Debug, Clone, PartialEq)]
124struct CurrentStatistics {
125    num_frames: u64,
126    false_decodes: u64,
127    total_iterations: u64,
128    start: Instant,
129    ldpc: CurrentCodeStatistics,
130    bch: Option<CurrentCodeStatistics>,
131}
132
133#[derive(Debug, Clone, PartialEq)]
134struct CurrentCodeStatistics {
135    bit_errors: u64,
136    frame_errors: u64,
137    correct_iterations: u64,
138}
139
140/// BER test statistics.
141///
142/// This structure contains the statistics for a single Eb/N0 case in a BER
143/// test.
144#[derive(Debug, Clone, PartialEq)]
145pub struct Statistics {
146    /// Eb/N0 in dB units.
147    pub ebn0_db: f32,
148    /// Number of frames tested.
149    pub num_frames: u64,
150    /// Total number of iterations.
151    pub total_iterations: u64,
152    /// Number of frames falsely decoded.
153    ///
154    /// These are frames for which the decoder has converged to a valid
155    /// codeword, but the codeword is different from the transmitted codeword.
156    pub false_decodes: u64,
157    /// Average iterations per frame.
158    pub average_iterations: f64,
159    /// Elapsed time for this test case.
160    pub elapsed: Duration,
161    /// Throughput in Mbps (referred to information bits).
162    pub throughput_mbps: f64,
163    /// Statistics of the inner LDPC decoder.
164    pub ldpc: CodeStatistics,
165    /// Statistics of the combined inner LDPC decoder plus outer BCH decoder (if it exists).
166    pub bch: Option<CodeStatistics>,
167}
168
169/// BER test statistics for a particular code.
170///
171/// This is a substructure of [`Statistics`] that includes all the elements that
172/// depend on whether we are using only the inner LDPC decoder or LDPC plus
173/// BCH. `Statistics` has two instances of this structure: one for LDPC-only and
174/// another for LDPC plus BCH (which is only present when BCH is enabled).
175#[derive(Debug, Clone, PartialEq)]
176pub struct CodeStatistics {
177    /// Number of bit errors.
178    pub bit_errors: u64,
179    /// Number of frame errors.
180    pub frame_errors: u64,
181    /// Sum of iterations in correct frames.
182    pub correct_iterations: u64,
183    /// Bit error rate.
184    pub ber: f64,
185    /// Frame error rate.
186    pub fer: f64,
187    /// Average iterations per correct frame.
188    pub average_iterations_correct: f64,
189}
190
191/// Progress reporter.
192///
193/// A reporter can optionally be supplied to the BER test on contruction in
194/// order to receive periodic messages reporting the test progress.
195#[derive(Debug, Clone)]
196pub struct Reporter {
197    /// Sender element of a channel used to send the reports.
198    pub tx: Sender<Report>,
199    /// Reporting interval.
200    pub interval: Duration,
201}
202
203/// BER test progress report.
204///
205/// Progress reports are optionally sent out periodically by the BER test. These
206/// can be used to update a UI to show the progress.
207#[derive(Debug, Clone, PartialEq)]
208pub enum Report {
209    /// Statistics for the current Eb/N0 being tested.
210    ///
211    /// This is sent periodically, and also when the Eb/N0 is finished.
212    Statistics(Statistics),
213    /// The complete BER test has finished.
214    ///
215    /// This is sent when all the Eb/N0 cases have been done.
216    Finished,
217}
218
219macro_rules! report {
220    ($self:expr, $current_statistics:expr, $ebn0_db:expr, $final:expr) => {
221        if let Some(reporter) = $self.reporter.as_ref() {
222            let now = Instant::now();
223            if $final || $self.last_reported + reporter.interval < now {
224                reporter
225                    .tx
226                    .send(Report::Statistics(Statistics::from_current(
227                        &$current_statistics,
228                        $ebn0_db,
229                        $self.k,
230                    )))
231                    .unwrap();
232                $self.last_reported = now;
233            }
234        }
235    };
236}
237
238impl<Mod: Modulation, Dec: DecoderFactory> BerTest<Mod, Dec> {
239    /// Creates a new BER test.
240    ///
241    /// The parameters required to define the test are given in the
242    /// [`BerTestParameters`] struct.
243    ///
244    /// This function only defines the BER test. To run it it is necessary to
245    /// call the [`BerTest::run`] method.
246    pub fn new(parameters: BerTestParameters<'_, Dec>) -> Result<BerTest<Mod, Dec>, Error> {
247        let k = parameters.h.num_cols() - parameters.h.num_rows();
248        let n_cw = parameters.h.num_cols();
249        let puncturer = parameters.puncturing_pattern.map(Puncturer::new);
250        let interleaver = parameters
251            .interleaving_columns
252            .map(|n| Interleaver::new(n.unsigned_abs(), n < 0));
253        let puncturer_rate = if let Some(p) = puncturer.as_ref() {
254            p.rate()
255        } else {
256            1.0
257        };
258        let n = (n_cw as f64 / puncturer_rate).round() as usize;
259        let rate = k as f64 / n as f64;
260        Ok(BerTest {
261            decoder_implementation: parameters.decoder_implementation,
262            num_workers: parameters.num_workers,
263            k,
264            n,
265            n_cw,
266            rate,
267            encoder: Encoder::from_h(&parameters.h)?,
268            h: parameters.h,
269            puncturer,
270            interleaver,
271            modulator: Mod::Modulator::default(),
272            ebn0s_db: parameters.ebn0s_db.to_owned(),
273            statistics: Vec::with_capacity(parameters.ebn0s_db.len()),
274            bch_max_errors: parameters.bch_max_errors,
275            max_iterations: parameters.max_iterations,
276            max_frame_errors: parameters.max_frame_errors,
277            min_run_time: parameters.min_run_time.unwrap_or(Duration::ZERO),
278            max_run_time: parameters.max_run_time.unwrap_or(Duration::MAX),
279            reporter: parameters.reporter,
280            last_reported: Instant::now(),
281        })
282    }
283
284    /// Runs the BER test.
285    ///
286    /// This function runs the BER test until completion. It returns a list of
287    /// statistics for each Eb/N0, or an error.
288    pub fn run(mut self) -> Result<Vec<Statistics>, Box<dyn std::error::Error>> {
289        let ret = self.do_run();
290        if let Some(reporter) = self.reporter.as_ref() {
291            reporter.tx.send(Report::Finished).unwrap();
292        }
293        ret?;
294        Ok(self.statistics)
295    }
296
297    fn do_run(&mut self) -> Result<(), Box<dyn std::error::Error>> {
298        self.last_reported = Instant::now();
299        for &ebn0_db in &self.ebn0s_db {
300            let ebn0 = 10.0_f64.powf(0.1 * f64::from(ebn0_db));
301            let esn0 = self.rate * Mod::BITS_PER_SYMBOL * ebn0;
302            let noise_sigma = (0.5 / esn0).sqrt();
303            let (results_tx, results_rx) = mpsc::channel();
304            let workers = std::iter::repeat_with(|| {
305                let (mut worker, terminate_tx) = self.make_worker(noise_sigma, results_tx.clone());
306                let handle = std::thread::spawn(move || worker.work());
307                (handle, terminate_tx)
308            })
309            .take(self.num_workers)
310            .collect::<Vec<_>>();
311
312            let mut current_statistics = CurrentStatistics::new(self.bch_max_errors > 0);
313            while !current_statistics.run_finished(
314                self.max_frame_errors,
315                self.min_run_time,
316                self.max_run_time,
317            ) {
318                match results_rx.recv().unwrap() {
319                    Ok(result) => {
320                        current_statistics.ldpc.bit_errors += result.bit_errors;
321                        current_statistics.ldpc.frame_errors += u64::from(result.frame_error);
322                        current_statistics.false_decodes += u64::from(result.false_decode);
323                        current_statistics.total_iterations += result.iterations;
324                        if !result.frame_error {
325                            current_statistics.ldpc.correct_iterations += result.iterations;
326                        }
327                        current_statistics.num_frames += 1;
328                        if let Some(bch) = &mut current_statistics.bch {
329                            if result.bit_errors > self.bch_max_errors {
330                                // BCH cannot decode codeword
331                                bch.bit_errors += result.bit_errors;
332                                bch.frame_errors += 1;
333                            } else {
334                                // BCH can decode codeword
335                                bch.correct_iterations += result.iterations;
336                            }
337                        }
338                    }
339                    Err(()) => break,
340                }
341                report!(self, current_statistics, ebn0_db, false);
342            }
343            report!(self, current_statistics, ebn0_db, true);
344
345            for (_, terminate_tx) in workers.iter() {
346                // we don't care if this fails because the worker has terminated
347                // and dropped the channel.
348                let _ = terminate_tx.send(());
349            }
350
351            let mut join_error = None;
352            for (handle, _) in workers.into_iter() {
353                if let Err(e) = handle.join().unwrap() {
354                    join_error = Some(e);
355                }
356            }
357            if let Some(e) = join_error {
358                return Err(e);
359            }
360
361            self.statistics.push(Statistics::from_current(
362                &current_statistics,
363                ebn0_db,
364                self.k,
365            ));
366        }
367        Ok(())
368    }
369
370    fn make_worker(
371        &self,
372        noise_sigma: f64,
373        results_tx: Sender<WorkerResult>,
374    ) -> (Worker<Mod>, SyncSender<()>) {
375        let (terminate_tx, terminate_rx) = mpsc::sync_channel(1);
376        (
377            Worker {
378                terminate_rx,
379                results_tx,
380                k: self.k,
381                encoder: self.encoder.clone(),
382                puncturer: self.puncturer.clone(),
383                interleaver: self.interleaver.clone(),
384                modulator: self.modulator.clone(),
385                channel: AwgnChannel::new(noise_sigma),
386                demodulator: Mod::Demodulator::from_noise_sigma(noise_sigma),
387                decoder: self.decoder_implementation.build_decoder(self.h.clone()),
388                max_iterations: self.max_iterations,
389            },
390            terminate_tx,
391        )
392    }
393}
394
395impl<Mod: Modulation, Dec: DecoderFactory> Ber for BerTest<Mod, Dec> {
396    fn run(self: Box<Self>) -> Result<Vec<Statistics>, Box<dyn std::error::Error>> {
397        BerTest::run(*self)
398    }
399
400    fn n(&self) -> usize {
401        self.n
402    }
403
404    fn n_cw(&self) -> usize {
405        self.n_cw
406    }
407
408    fn k(&self) -> usize {
409        self.k
410    }
411
412    fn rate(&self) -> f64 {
413        self.rate
414    }
415}
416
417impl<Mod: Modulation> Worker<Mod> {
418    fn work(&mut self) -> Result<(), Box<dyn std::error::Error + Send + Sync + 'static>> {
419        let mut rng = rand::rng();
420        loop {
421            match self.terminate_rx.try_recv() {
422                Ok(()) => return Ok(()),
423                Err(TryRecvError::Disconnected) => panic!(),
424                Err(TryRecvError::Empty) => (),
425            };
426            let result = self.simulate(&mut rng);
427            let to_send = match result.as_ref() {
428                Ok(r) => Ok(r.clone()),
429                Err(_) => Err(()),
430            };
431            self.results_tx.send(to_send).unwrap();
432            result?;
433        }
434    }
435
436    fn simulate<R: Rng>(
437        &mut self,
438        rng: &mut R,
439    ) -> Result<WorkerResultOk, Box<dyn std::error::Error + Send + Sync + 'static>> {
440        let message = Self::random_message(rng, self.k);
441        let codeword = self.encoder.encode(&Self::gf2_array(&message));
442        let transmitted = match self.puncturer.as_ref() {
443            Some(p) => p.puncture(&codeword)?,
444            None => codeword,
445        };
446        let transmitted = match self.interleaver.as_ref() {
447            Some(i) => i.interleave(&transmitted),
448            None => transmitted,
449        };
450        let mut symbols = self.modulator.modulate(&transmitted);
451        self.channel.add_noise(rng, &mut symbols);
452        let llrs_demod = self.demodulator.demodulate(&symbols);
453        let llrs_decoder = match self.interleaver.as_ref() {
454            Some(i) => i.deinterleave(&llrs_demod),
455            None => llrs_demod,
456        };
457        let llrs_decoder = match self.puncturer.as_ref() {
458            Some(p) => p.depuncture(&llrs_decoder)?,
459            None => llrs_decoder,
460        };
461
462        let (decoded, iterations, success) =
463            match self.decoder.decode(&llrs_decoder, self.max_iterations) {
464                Ok(output) => (output.codeword, output.iterations, true),
465                Err(output) => (output.codeword, output.iterations, false),
466            };
467        // Count only bit errors in the systematic part of the codeword
468        let bit_errors = message
469            .iter()
470            .zip(decoded.iter())
471            .filter(|(a, b)| a != b)
472            .count() as u64;
473        let frame_error = bit_errors > 0;
474        let false_decode = frame_error && success;
475        Ok(WorkerResultOk {
476            bit_errors,
477            frame_error,
478            false_decode,
479            iterations: iterations as u64,
480        })
481    }
482
483    fn random_message<R: Rng>(rng: &mut R, size: usize) -> Vec<u8> {
484        rng.sample_iter(StandardUniform)
485            .map(<u8 as From<bool>>::from)
486            .take(size)
487            .collect()
488    }
489
490    fn gf2_array(bits: &[u8]) -> Array1<GF2> {
491        Array1::from_iter(
492            bits.iter()
493                .map(|&b| if b == 1 { GF2::one() } else { GF2::zero() }),
494        )
495    }
496}
497
498impl CurrentStatistics {
499    fn new(has_bch: bool) -> CurrentStatistics {
500        CurrentStatistics {
501            num_frames: 0,
502            false_decodes: 0,
503            total_iterations: 0,
504            start: Instant::now(),
505            ldpc: CurrentCodeStatistics::new(),
506            bch: if has_bch {
507                Some(CurrentCodeStatistics::new())
508            } else {
509                None
510            },
511        }
512    }
513
514    fn errors_for_termination(&self) -> u64 {
515        if let Some(bch) = &self.bch {
516            bch.frame_errors
517        } else {
518            self.ldpc.frame_errors
519        }
520    }
521
522    fn run_finished(
523        &self,
524        max_frame_errors: u64,
525        min_run_time: Duration,
526        max_run_time: Duration,
527    ) -> bool {
528        let elapsed = self.start.elapsed();
529        (self.errors_for_termination() >= max_frame_errors && elapsed >= min_run_time)
530            || elapsed >= max_run_time
531    }
532}
533
534impl CurrentCodeStatistics {
535    fn new() -> CurrentCodeStatistics {
536        CurrentCodeStatistics {
537            bit_errors: 0,
538            frame_errors: 0,
539            correct_iterations: 0,
540        }
541    }
542}
543
544impl Default for CurrentCodeStatistics {
545    fn default() -> CurrentCodeStatistics {
546        CurrentCodeStatistics::new()
547    }
548}
549
550impl CodeStatistics {
551    fn from_current(stats: &CurrentCodeStatistics, num_frames: u64, k: usize) -> CodeStatistics {
552        CodeStatistics {
553            bit_errors: stats.bit_errors,
554            frame_errors: stats.frame_errors,
555            correct_iterations: stats.correct_iterations,
556            ber: stats.bit_errors as f64 / (k as f64 * num_frames as f64),
557            fer: stats.frame_errors as f64 / num_frames as f64,
558            average_iterations_correct: stats.correct_iterations as f64
559                / (num_frames - stats.frame_errors) as f64,
560        }
561    }
562}
563
564impl Statistics {
565    fn from_current(stats: &CurrentStatistics, ebn0_db: f32, k: usize) -> Statistics {
566        let elapsed = Instant::now() - stats.start;
567        Statistics {
568            ebn0_db,
569            num_frames: stats.num_frames,
570            false_decodes: stats.false_decodes,
571            total_iterations: stats.total_iterations,
572            average_iterations: stats.total_iterations as f64 / stats.num_frames as f64,
573            elapsed,
574            throughput_mbps: 1e-6 * (k as f64 * stats.num_frames as f64) / elapsed.as_secs_f64(),
575            ldpc: CodeStatistics::from_current(&stats.ldpc, stats.num_frames, k),
576            bch: stats
577                .bch
578                .as_ref()
579                .map(|b| CodeStatistics::from_current(b, stats.num_frames, k)),
580        }
581    }
582}
ldpc_toolbox/simulation/ber.rs

ldpc_toolbox/simulation/
ber.rs
