new test program

diff --git a/test/test.c b/test/test.c
new file mode 100644 (file)
index 0000000..6e9b5a3
--- /dev/null
+++ b/test/test.c
@@ -0,0 +1,245 @@
+/* -*- C -*- */
+
+/*
+
+Copyright (C) 2025 Laurent Mazet
+
+*/
+
+/* Program test for Viterbi decoding */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "trellis.h"
+#include "encoder.h"
+#include "create_trellis.h"
+
+/* Pseudo-random bit sequence */
+void bit_gene (bit * a, int n)
+{
+    for (int i = 0; i < n; i++) {
+        a[i] = rand () & 1;
+    }
+}
+
+#if !(_SVID_SOURCE || _XOPEN_SOURCE)
+double drand48 (void) {
+    return rand() / (RAND_MAX + 1.0);
+}
+#endif
+
+/* Gaussian random variable */
+double randn (void)
+{
+    static int iset = 0;
+    static double gset = 0;
+
+    if (iset) {
+        iset = 0;
+        return gset;
+    }
+
+    double v1, v2, rsq;
+    do {
+      v1 = 2.0*drand48()-1.0;
+      v2 = 2.0*drand48()-1.0;
+      rsq = v1*v1+v2*v2;
+    } while ((rsq >= 1) || (rsq == 0));
+    double fac = sqrt (-2 * log (rsq) / rsq);
+    gset = v1*fac;
+    iset = 1;
+
+    return v2 * fac;
+}
+
+int main (int argc, char *argv[])
+{
+
+    /* Coders definition */
+    int nb_unprotected_bits = 0;
+    int nb_coded_inputs = 1;
+    int nb_coded_outputs = 2;
+    int nb_registers = 2;
+
+    /* Data paramaters */
+    int nb_data_bits = 128 * 8;
+    int nb_packets = 1024 * 10;
+
+    /* Argument processing */
+    if ((argc > 1) && (argv[1][0] == '-') && (argv[1][1] == 'h')) {
+        printf ("prog_test [-h] [#unprotected #inputs #outputs #registers"
+                " coder_1 coder_2...]\n");
+        return 0;
+    }
+    if (argc > 1) {
+        nb_unprotected_bits = atoi (argv[1]);
+    }
+    if (argc > 2) {
+        nb_coded_inputs = atoi (argv[2]);
+    }
+    if (argc > 3) {
+        nb_coded_outputs = atoi (argv[3]);
+    }
+    if (argc > 4) {
+        nb_registers = atoi (argv[4]);
+    }
+
+    int *coders = (int *) calloc (nb_coded_outputs, sizeof (int));
+    if (argc == 1) {
+        coders[0] = 5;
+        coders[1] = 7;
+    } else if (argc == 5 + nb_coded_outputs) {
+        for (int i = 0; i < nb_coded_outputs; i++) {
+            coders[i] = (int) strtol (argv[5 + i], NULL, 8);
+        }
+    } else {
+        printf ("Wrong number of arguments\n");
+        return 1;
+    }
+
+    /* Encoder descriptions */
+    printf ("#S/#I/#O/#R: (%d,%d,%d,%d) [", nb_unprotected_bits,
+            nb_coded_inputs, nb_coded_outputs, nb_registers);
+    for (int i = 0; i < nb_coded_outputs; i++) {
+        printf (" %o", coders[i]);
+    }
+    printf (" ]\n");
+
+    /* Coders definition */
+    int nb_inputs = nb_unprotected_bits + nb_coded_inputs;
+    int nb_outputs = nb_unprotected_bits + nb_coded_outputs;
+    conv_encoder *enc =
+        alloc_conv_encoder (nb_unprotected_bits, nb_coded_inputs, nb_registers,
+                            nb_coded_outputs, coders);
+
+    /* Data paramaters */
+    int nb_uncoded_bits = nb_data_bits + nb_registers;
+    int nb_coded_bits = nb_uncoded_bits * nb_outputs;
+    int nb_metrics = 1 << nb_outputs;
+
+    /* Viterbi parameters for a closed trellis */
+    int deepness = nb_registers / nb_inputs;
+    int window_dec_size = nb_uncoded_bits / nb_inputs - deepness;
+    int history_size = deepness + window_dec_size;
+
+    /* Trellis creation */
+    trellis *tr = create_trellis_from_conv_encoder (enc, history_size);
+
+    /* Info */
+    printf ("Send %d packets of %d bits\n", nb_packets, nb_data_bits);
+
+    for (double EbN0 = -15; EbN0 < 50; EbN0 += 0.5) {
+        double sigma = 1. / sqrt (2.) * pow (10., -EbN0 / 20.);
+    
+        /* Loop for benchmark */
+        int nb_total_bit_errors = 0;
+        int nb_packet_errors = 0;
+        int nb_packet_send = 0;
+        for (int l = 0; l < nb_packets; l++) {
+    
+            /* Some malloc */
+            bit *x = (bit *) calloc (nb_uncoded_bits, sizeof (bit));
+            bit *y = (bit *) calloc (nb_coded_bits, sizeof (bit));
+            metric_t *metrics = (metric_t *) calloc (nb_metrics, sizeof (metric_t));
+            bit *xdecod = (bit *) calloc (nb_uncoded_bits, sizeof (bit));
+    
+            /* Bit generation */
+            bit_gene (x, nb_data_bits);
+            for (int i = nb_data_bits; i < nb_uncoded_bits; i++) {
+                x[i] = 0;
+            }
+    
+            /* Bit encoding */
+            conv_encoding (y, x, nb_uncoded_bits, enc);
+    
+            /* Viterbi */
+            for (int i = 0, ii = 0; i < nb_uncoded_bits; i += nb_inputs) {
+                int symbol = 0;
+    
+                /* Generate symbol */
+                for (int k = 0; k < nb_outputs; k++) {
+                    symbol = (symbol << 1) | y[ii++];
+                }
+    
+                /* Compute distance of received signal to all constellation symbols */
+                for (int k = 0; k < nb_metrics; k++) {
+                    int bits = symbol ^ k;
+    
+                    int m = 0;
+                    for (int l = 0; (l < nb_outputs) && (bits); l++, bits >>= 1) {
+                        if (bits & 1) {
+                            m++;
+                        }
+                    }
+                    metrics[k] = m + randn() * sigma;
+                }
+                /* acs */
+                tr->acs (tr, metrics);
+    
+                /* Normalize node metrics */
+                /* normalize_node_metrics_trellis (tr, 0); */
+            }
+    
+            /* Decoded bits */
+            int *tmp = (int *) calloc (window_dec_size, sizeof (int));
+            int nb_dec_symb = backtrace_trellis (tr, 0, deepness, window_dec_size, tmp);
+    
+            for (int j = 0; j < nb_dec_symb; j++) {
+                for (int k = 0; k < nb_inputs; k++) {
+                    xdecod[j * nb_inputs + k] =
+                        ((tmp[j] & (1 << (nb_inputs - k - 1))) != 0) ? 1 : 0;
+                }
+            }
+            free (tmp);
+    
+            /* Check error */
+            int nb_bit_errors = 0;
+            for (int i = 0; i < nb_data_bits; i++) {
+                if (x[i] != xdecod[i]) {
+                    nb_bit_errors++;
+                }
+            }
+    
+            /* Free memory blocks */
+            free (x);
+            free (y);
+            free (metrics);
+            free (xdecod);
+    
+            reset_trellis (tr);
+    
+            if (nb_bit_errors) {
+                //printf ("Nb bit errors: %d\n", nb_bit_errors);
+                nb_packet_errors++;
+            }
+            nb_total_bit_errors += nb_bit_errors;
+            nb_packet_send++;
+            if (nb_packet_errors == 10) {
+                break;
+            }
+        }
+    
+        double ber = (double)nb_total_bit_errors / (nb_packet_send * nb_data_bits);
+        double per = (double)nb_packet_errors / nb_packet_send;
+        char alert = (nb_packet_send == nb_packets) ? '!' : '\0';
+        if (nb_packet_errors != nb_packet_send) {
+            printf ("EbN0: %5.1f BER: %.1e PER: %.1e%c\n", EbN0, ber, per, alert);
+        }
+    
+        if (per < 1e-2) {
+            break;
+        }
+    }
+
+    free_trellis (tr);
+    free_conv_encoder (enc);
+    free (coders);
+
+    printf ("bye\n");
+
+    return 0;
+}
+
+/* vim: set ts=4 sw=4 si et: */