Viterbi

Bakgrund

Viterbis algoritm – Hur vi hittar rätt signal i bruset

Inom amatörradion och digital signalbehandling stöter vi ofta på problemet med brus och fading. När vi skickar digital data över etern, till exempel med modernare digitala trafiksätt eller vid satellitkommunikation, använder vi ofta något som kallas för faltningskodning (convolutional coding). Det är en metod där varje sänd bit beror på både den nuvarande biten och ett antal tidigare bitar. Detta skapar ett mönster, ett beroende, som gör signalen robust.

Men hur avkodar vi detta på mottagarsidan när atmosfären har ställt till det och lagt till en massa brus? Det är här Viterbis algoritm kommer in i bilden.

Kort sagt är Viterbialgoritmen en metod för att hitta den mest troliga sekvensen av dolda tillstånd (det vi faktiskt sände) baserat på en serie observerade händelser (den brusiga signalen vi tog emot). Istället för att gissa bit för bit, tittar algoritmen på hela sammanhanget.

Trellis-diagrammet: Vägarna genom historien

För att förstå hur algoritmen jobbar kan vi tänka oss ett nätverk av vägar, ett så kallat trellis-diagram.

• Tillstånd (States): Radiosändarens "minne" kan vid varje given tidpunkt befinna sig i ett visst antal tillstånd (beroende på de senaste sända bitarna).

• Övergångar: När en ny bit ska sändas hoppar sändaren till ett nytt tillstånd och skickar ut en specifik kombination av bitar i luften.

Hur funkar det?

När mottagaren lyssnar får den in en sekvens som kanske är halvt dränkt i brus. En nolla kanske ser ut som en svag etta, och så vidare. Viterbialgoritmen börjar nu räkna på alla möjliga vägar genom detta nätverk av tillstånd över tid. Så här jobbar algoritmen steg för steg:

1. Beräkna avstånd (Metric): För varje steg i tiden tittar algoritmen på de bitar som faktiskt togs emot och jämför dem med vad som borde ha sänts för varje tänkbar väg i diagrammet. Man mäter "avståndet" (ofta Hamming-avstånd för ettor och nollor, eller Euklidiskt avstånd om man kör "soft decision" där man mäter signalstyrkan mer exakt). Ju mindre avstånd, desto mer troligt är det att den vägen är den rätta.

1. Spara den bästa vägen (Add-Compare-Select): Till varje nytt tillstånd i diagrammet leder det oftast två olika vägar från det föregående steget. Algoritmen adderar det nya avståndet till det gamla, jämför de två alternativen, och väljer den väg som har lägst totalt avstånd (fel). Den andra vägen slängs bort för alltid. Det här är algoritmens stora styrka: den rensar hela tiden bort återvändsgränder så att beräkningsbördan inte exploderar.

1. Lås historiken (Traceback): När hela sekvensen (eller ett tillräckligt långt block) har tagits emot, tittar man på vilket slutgiltigt tillstånd som har det absolut lägsta totala felvärdet. Sedan backar man bakåt genom de sparade "bästa vägarna" till början.

Resultatet blir den absolut mest sannolika sekvensen av bitar som sändaren faktiskt skickade ut, trots att det var fullt av bitfel under resans gång.

Implementation

Modern C

/*
 * Anders Sikvall 2026 Version 1.0
 * <anders@sikvall.se>
 * 
 * Implementation i C av en standard Viterbi-algoritm som används bl.a. av 
 * Qualcomm och NASA och är välpublicerad.
 *
 * PUBLIC DOMAIN FOR EDUCATIONAL PURPOSES
 *
 * Kompilera med gcc -o viterbi viterbi.c
 *
 * ANVÄNDNING
 *         viterbi -e <indatastream >utdatastream   Kodar <indatastream till >utdatastream
 *         viterbi -d <indatastream >utdatastream   Avkodar <indatastream till >utdatastream
 *         viterbi -t                               Utför en intern test på 64 randomiserade oktetter
 *         viterbi -h                               Visar en kort hjälptext
 *         viterbi --error                          Samma sak som -t men introducerar 3 flippade bitpar
 */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <stdint.h>
#include <time.h>

/* Constraint length (K=7) innebär att skiftregistret minns de senaste 7 bitarna. */
#define K 7

/* Antalet tillstånd i trellis-diagrammet bestäms av de senaste K-1 bitarna.
   2^(7-1) = 2^6 = 64 unika tillstånd totalt. */
#define NUM_STATES (1 << (K - 1))

/* Standard NASA/Qualcomm-polynom för Rate 1/2 konvolutionskodning.
   POLY_G0 motsvarar oktalt 151 (1001111 binärt).
   POLY_G1 motsvarar oktalt 113 (1101101 binärt).
   Dessa maskar används för att bestämma vilka bitar i tillståndet som XOR-sammanfogas. */
const uint8_t POLY_G0 = 0x4F;
const uint8_t POLY_G1 = 0x6D;

/* Hårdvaruoptimerad paritetsberäkning.
   Räknar antalet satta bitar (1:or) i ett tal och returnerar 1 om antalet är udda,
   eller 0 om antalet är jämnt (motsvarar en kedja av XOR-grindar). */
static inline uint8_t parity(uint8_t x) {
    return __builtin_parity(x) & 1;
}

// ============================================================================
// KODARE (-e)
// ============================================================================
void encode(FILE *in, FILE *out) {
    /* encoder_state lagrar historiken av de senaste K-1 (6) bitarna. */
    uint8_t encoder_state = 0;

    /* out_byte används som en tillfällig buffer för att samla ihop bitar till en hel byte. */
    uint8_t out_byte = 0;

    /* bit_count håller koll på hur många bitar som för tillfället ligger i utmatningsbufferten. */
    int bit_count = 0;
    int ch;

    /* Loopa igenom infilen tecken för tecken (byte för byte). */
    while ((ch = fgetc(in)) != EOF) {

        /* Gå igenom varje enskild bit i den inlästa byten, från MSB (bit 7) till LSB (bit 0). */
        for (int i = 7; i >= 0; i--) {

            /* Isolera biten på position i. */
            uint8_t bit = (ch >> i) & 1;

            /* Skapa det nuvarande ordet genom att skifta upp det gamla tillståndet till vänster
               och lägga till den nya inkommande biten som den lägsta biten (LSB). */
            uint8_t current_word = (encoder_state << 1) | bit;

            /* Beräkna de två paritetsbitarna (out0 och out1) genom att applicera polynommaskerna
               och köra XOR (paritet) på resultatet. */
            uint8_t out0 = parity(current_word & POLY_G0);
            uint8_t out1 = parity(current_word & POLY_G1);

            /* Knuffa in den första paritetsbiten (out0) i utmatningsbyten. */
            out_byte = (out_byte << 1) | out0;
            bit_count++;

            /* Om vi har fyllt en hel byte (8 bitar), skriv ut den till utmatningsströmmen. */
            if (bit_count == 8) {
                fputc(out_byte, out);
                out_byte = 0;
                bit_count = 0;
            }

            /* Knuffa in den andra paritetsbiten (out1) i utmatningsbyten. */
            out_byte = (out_byte << 1) | out1;
            bit_count++;

            /* Kontrollera återigen om utmatningsbyten blev full efter den andra biten. */
            if (bit_count == 8) {
                fputc(out_byte, out);
                out_byte = 0;
                bit_count = 0;
            }

            /* Uppdatera kodarens interna tillstånd inför nästa bit.
               Vi sparar endast de lägsta K-1 bitarna med hjälp av en bitmask (63 = 0x3F). */
            encoder_state = current_word & (NUM_STATES - 1);
        }
    }

    /* Om det finns kvarvarande bitar i bufferten när filen slutar,
       skiftar vi dem till vänster så att de ligger rätt, och skriver ut den sista ofullständiga byten. */
    if (bit_count > 0) {
        out_byte <<= (8 - bit_count);
        fputc(out_byte, out);
    }
}

// ============================================================================
// AVKODARE (-d)
// ============================================================================
void decode(FILE *in, FILE *out) {
    /* Startkapacitet för historikmatrisen (antalet bit-steg). Den växer linjärt vid behov. */
    size_t history_capacity = 1024;
    size_t step = 0;

    /* Allokera historikmatrisen dynamiskt. Det är en 2D-array av storleken [kapacitet][64].
       Den sparar överlevnadsbesluten (0 eller 1) för varje tillstånd vid varje givet tidssteg. */
    uint8_t (*history)[NUM_STATES] = malloc(history_capacity * sizeof(*history));
    if (!history) {
        return;
    }

    /* path_metrics lagrar den ackumulerade felkostnaden (Hammingavståndet) för de 64 nuvarande tillstånden. */
    uint32_t path_metrics[NUM_STATES];

    /* old_metrics lagrar en oförändrad kopia av kostnaderna från föregående tidssteg. */
    uint32_t old_metrics[NUM_STATES];

    /* Initiera alla tillstånd till en mycket hög felkostnad ("oändlighet").
       Undantaget är tillstånd 0, vilket är det kända och garanterade starttillståndet (kostnad 0). */
    for (int i = 1; i < NUM_STATES; i++) {
        path_metrics[i] = 1e6;
    }
    path_metrics[0] = 0;

    int ch;

    // ------------------------------------------------------------------------
    // STEG 1: Framåtpass (Forward Pass / Add-Compare-Select)
    // ------------------------------------------------------------------------
    /* Läs den kodade strömmen byte för byte. */
    while ((ch = fgetc(in)) != EOF) {

        /* Varje inläst byte innehåller 4 bitpar (totalt 8 kodade bitar), vilket motsvarar 4 tidssteg. */
        for (int bit_pair = 3; bit_pair >= 0; bit_pair--) {

            /* Extrahera de två bitarna (b0 och b1) som hör till det aktuella steget. */
            uint8_t b0 = (ch >> (bit_pair * 2 + 1)) & 1;
            uint8_t b1 = (ch >> (bit_pair * 2)) & 1;

            /* Om vår historikmatris är full, fördubblar vi dess storlek dynamiskt med realloc.
               Detta gör att programmet kan hantera obegränsat stora filer och pipes. */
            if (step >= history_capacity) {
                history_capacity *= 2;
                uint8_t (*new_history)[NUM_STATES] = realloc(history, history_capacity * sizeof(*history));
                if (!new_history) {
                    free(history);
                    return;
                }
                history = new_history;
            }

            /* Kopiera över nuvarande path_metrics till old_metrics inför beräkningen av nästa steg. */
            memcpy(old_metrics, path_metrics, sizeof(path_metrics));

            /* Loopa igenom alla 64 möjliga måltillstånd (nästa tillstånd) i trellis-diagrammet. */
            for (int state = 0; state < NUM_STATES; state++) {

                /* Eftersom vi skiftar till vänster, kan vi räkna ut vilka två föregående tillstånd
                   (prev0 och prev1) som kan ha skiftat in i det nuvarande tillståndet. */
                int prev0 = state >> 1;
                int prev1 = prev0 | (NUM_STATES >> 1);

                /* Den inbit som krävdes för att nå 'state' är identisk med nuvarande tillstånds LSB. */
                uint8_t input_bit = state & 1;

                /* --- Beräkna förväntade bitar och Hammingavstånd för VÄG 0 (från prev0) --- */
                uint8_t word0 = (prev0 << 1) | input_bit;
                uint8_t exp0_0 = parity(word0 & POLY_G0);
                uint8_t exp0_1 = parity(word0 & POLY_G1);
                /* Kostnaden ökar med 1 för varje bit som inte matchar de faktiskt mottagna bitarna. */
                uint32_t cost0 = (exp0_0 != b0) + (exp0_1 != b1);

                /* --- Beräkna förväntade bitar och Hammingavstånd för VÄG 1 (från prev1) --- */
                uint8_t word1 = (prev1 << 1) | input_bit;
                uint8_t exp1_0 = parity(word1 & POLY_G0);
                uint8_t exp1_1 = parity(word1 & POLY_G1);
                uint32_t cost1 = (exp1_0 != b0) + (exp1_1 != b1);

                /* Addera det gamla ackumulerade felet med det nya övergångsfelet för båda vägarna. */
                uint32_t m0 = old_metrics[prev0] + cost0;
                uint32_t m1 = old_metrics[prev1] + cost1;

                /* (Compare & Select): Välj den väg som har lägst totalt fel.
                   Spara det nya lägsta felet i path_metrics, och dokumentera valet (0 eller 1)
                   i historikmatrisen för att kunna återskapa vägen baklänges senare. */
                if (m0 <= m1) {
                    path_metrics[state] = m0;
                    history[step][state] = 0;
                } else {
                    path_metrics[state] = m1;
                    history[step][state] = 1;
                }
            }
            /* Öka räknaren för antal avklarade bit-steg. */
            step++;
        }
    }

    /* Om vi inte lyckades läsa några bitar alls, städa upp och avsluta. */
    if (step == 0) {
        free(history);
        return;
    }

    // ------------------------------------------------------------------------
    // STEG 2: Hitta bästa sluttillstånd
    // ------------------------------------------------------------------------
    /* Sök igenom alla 64 tillstånd vid slutet av strömmen för att hitta det tillstånd
       som har samlat på sig minst antal fel (lägst metric). */
    int curr_state = 0;
    uint32_t min_m = path_metrics[0];
    for (int i = 1; i < NUM_STATES; i++) {
        if (path_metrics[i] < min_m) {
            min_m = path_metrics[i];
            curr_state = i;
        }
    }

    // ------------------------------------------------------------------------
    // STEG 3: Bakåtspårning (Traceback)
    // ------------------------------------------------------------------------
    /* Allokera en temporär array för att spara de avkodade bitarna.
       Eftersom vi vandrar baklänges kommer vi att fylla denna array bakifrån. */
    uint8_t *decoded_bits = malloc(step);
    if (!decoded_bits) {
        free(history);
        return;
    }

    /* Vandra baklänges genom historiken från det sista tidssteget ner till noll. */
    for (int s = (int)step - 1; s >= 0; s--) {

        /* Den ursprungliga inbiten är alltid LSB (bit 0) i det tillstånd vi står i. */
        decoded_bits[s] = curr_state & 1;

        /* Hämta beslutet (0 eller 1) som fattades för detta tillstånd i tidssteg s. */
        int decision = history[s][curr_state];
        int prev0 = curr_state >> 1;

        /* Rekonstruera vilket tillstånd vi kom ifrån baserat på det sparade beslutet. */
        if (decision == 0) {
            curr_state = prev0;
        } else {
            curr_state = prev0 | (NUM_STATES >> 1);
        }
    }

    // ------------------------------------------------------------------------
    // STEG 4: Återpacka till bytes och skriv ut
    // ------------------------------------------------------------------------
    /* Bitarna ligger nu i korrekt kronologisk ordning i decoded_bits-arrayen.
       Vi loopar igenom dem och packar ihop dem 8 och 8 till vanliga bytes igen. */
    uint8_t current_byte = 0;
    for (size_t i = 0; i < step; i++) {
        current_byte = (current_byte << 1) | decoded_bits[i];

        /* Varje gång vi har samlat ihop 8 bitar skriver vi ut byten till utfilen/stdout. */
        if ((i + 1) % 8 == 0) {
            fputc(current_byte, out);
            current_byte = 0;
        }
    }

    /* Frigör det dynamiskt allokerade minnet. */
    free(decoded_bits);
    free(history);
}

// ============================================================================
// INTERNT TEST (-t OCH --error)
// ============================================================================
void run_test(bool inject_errors) {
    if (inject_errors) {
        printf("[*] Startar internt test med FELKORRIGERING (64 bytes slumpdata + 3 bitpar flippade)...\n");
    } else {
        printf("[*] Startar internt test (64 bytes randomiserat data)...\n");
    }

    /* Skapa en array med 64 bytes slumpmässigt data. */
    uint8_t original_data[64];
    srand((unsigned int)time(NULL));
    for (int i = 0; i < 64; i++) {
        original_data[i] = rand() % 256;
    }

    /* tmpfile() skapar temporära filer som körs i RAM/hårddisk och raderas automatiskt vid stängning. */
    FILE *orig_stream = tmpfile();
    FILE *encoded_stream = tmpfile();
    FILE *corrupted_stream = tmpfile();
    FILE *decoded_stream = tmpfile();

    if (!orig_stream || !encoded_stream || !corrupted_stream || !decoded_stream) {
        fprintf(stderr, "Kunde inte skapa temporära filer.\n");
        return;
    }

    /* Skriv originaldata till startströmmen och spola tillbaka pekaren till början. */
    fwrite(original_data, 1, 64, orig_stream);
    rewind(orig_stream);

    /* 1. Kör kodningen (64 bytes rådata blir exakt 128 bytes packad bitström). */
    encode(orig_stream, encoded_stream);
    rewind(encoded_stream);

    /* 2. Hantering av felinjicering (--error). */
    if (inject_errors) {
        uint8_t encoded_data[128];

        /* Läs in hela den kodade strömmen till en lokal array. */
        if (fread(encoded_data, 1, 128, encoded_stream) == 128) {

            /* Välj ut 3 slumpmässiga bytes och flippa en slumpmässig bit i var och en av dem. */
            for (int i = 0; i < 3; i++) {
                int byte_index = rand() % 128;
                int bit_index = rand() % 8;

                /* XOR-operatör (^) flippar den valda biten (0 blir 1, 1 blir 0). */
                encoded_data[byte_index] ^= (1 << bit_index);
                printf("[i] Introducerade bitfel i kodad byte %d, bit %d\n", byte_index, bit_index);
            }
        }
        /* Skriv ut det saboterade datat till den korrupta strömmen. */
        fwrite(encoded_data, 1, 128, corrupted_stream);
        rewind(corrupted_stream);
    } else {
        /* Om inga fel ska testas, pekar vi bara vidare till den felfria strömmen direkt. */
        corrupted_stream = encoded_stream;
    }

    /* 3. Kör avkodaren på dataströmmen. */
    decode(corrupted_stream, decoded_stream);
    rewind(decoded_stream);

    /* 4. Läs tillbaka det färdigavkodade datat och verifiera resultatet. */
    uint8_t final_data[64];
    size_t read_bytes = fread(final_data, 1, 64, decoded_stream);

    bool success = true;
    if (read_bytes != 64) {
        printf("[FAIL] Avkodade fel antal bytes: %zu (förväntade 64)\n", read_bytes);
        success = false;
    } else {
        /* Jämför ursprunglig array med den avkodade arrayen byte för byte. */
        for (int i = 0; i < 64; i++) {
            if (original_data[i] != final_data[i]) {
                printf("[FAIL] Diff vid byte %d: Org=0x%02X, Avkodad=0x%02X\n", i, original_data[i], final_data[i]);
                success = false;
                break;
            }
        }
    }

    /* Stäng alla strömmar ordentligt. */
    fclose(orig_stream);
    if (inject_errors) {
        fclose(encoded_stream);
        fclose(corrupted_stream);
    } else {
        fclose(encoded_stream);
    }
    fclose(decoded_stream);

    /* Skriv ut slutrapport för testet. */
    if (success) {
        if (inject_errors) {
            printf("[SUCCESS] Viterbi-algoritmens felkorrigering lyckades! Alla bitfel lagades och datat matchar perfekt.\n");
        } else {
            printf("[SUCCESS] Testet lyckades! Allt data matchar perfekt.\n");
        }
    }
}

// ============================================================================
// HJÄLPTEXT OCH MAIN
// ============================================================================
void print_help(const char *prog_name) {
    printf("Användning: %s [option]\n", prog_name);
    printf("Optioner:\n");
    printf("  -e <infil >utfil   Kodar infil till utfil (Rate 1/2, K=7)\n");
    printf("  -d <infil >utfil   Avkodar infil till utfil med Viterbi\n");
    printf("  -h                 Visar denna hjälptext\n");
    printf("  -t                 Kör internt test med 64 bytes slumpdata\n");
    printf("  --error            Kör internt test men injicerar bitfel för att testa felkorrigeringen\n");
}

int main(int argc, char *argv[]) {
    /* Säkerställ att användaren har skickat med exakt ett kommandoradsargument. */
    if (argc != 2) {
        print_help(argv[0]);
        return EXIT_FAILURE;
    }

    /* Kontrollera vilket argument användaren skickade och anropa rätt funktion. */
    if (strcmp(argv[1], "-e") == 0) {
        encode(stdin, stdout);
    } else if (strcmp(argv[1], "-d") == 0) {
        decode(stdin, stdout);
    } else if (strcmp(argv[1], "-h") == 0) {
        print_help(argv[0]);
    } else if (strcmp(argv[1], "-t") == 0) {
        run_test(false);
    } else if (strcmp(argv[1], "--error") == 0) {
        run_test(true);
    } else {
        fprintf(stderr, "Fel: Okänd parameter '%s'\n", argv[1]);
        print_help(argv[0]);
        return EXIT_FAILURE;
    }

    return EXIT_SUCCESS;
}