Viterbi

Från Täpp-Anders
Hoppa till navigeringHoppa till sök

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
 */


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

#define K 7
#define NUM_STATES (1 << (K - 1))

// Standard NASA/Qualcomm-polynom (Rate 1/2, K=7)
const uint8_t POLY_G0 = 0x4F; // 1001111 i binärt
const uint8_t POLY_G1 = 0x6D; // 1101101 i binärt

// Hårdvaruoptimerad paritetsberäkning
static inline uint8_t parity(uint8_t x) {
    return __builtin_parity(x) & 1;
}

// ============================================================================
// KODARE (-e)
// ============================================================================
void encode(FILE *in, FILE *out) {
    uint8_t encoder_state = 0;
    uint8_t out_byte = 0;
    int bit_count = 0;
    int ch;

    while ((ch = fgetc(in)) != EOF) {
        // Processa varje bit från MSB till LSB
        for (int i = 7; i >= 0; i--) {
            uint8_t bit = (ch >> i) & 1;
            uint8_t current_word = (encoder_state << 1) | bit;

            // Beräkna de två paritetsbitarna
            uint8_t out0 = parity(current_word & POLY_G0);
            uint8_t out1 = parity(current_word & POLY_G1);

            // Packa bit 0
            out_byte = (out_byte << 1) | out0;
            bit_count++;
            if (bit_count == 8) {
                fputc(out_byte, out);
                out_byte = 0;
                bit_count = 0;
            }

            // Packa bit 1
            out_byte = (out_byte << 1) | out1;
            bit_count++;
            if (bit_count == 8) {
                fputc(out_byte, out);
                out_byte = 0;
                bit_count = 0;
            }

            // Uppdatera skiftregistret (behåll K-1 bitar)
            encoder_state = current_word & (NUM_STATES - 1);
        }
    }
    
    // Töm resterande bitar i sista byten vid behov
    if (bit_count > 0) {
        out_byte <<= (8 - bit_count);
        fputc(out_byte, out);
    }
}

// ============================================================================
// AVKODARE (-d)
// ============================================================================
void decode(FILE *in, FILE *out) {
    size_t history_capacity = 1024;
    size_t step = 0;

    uint8_t (*history)[NUM_STATES] = malloc(history_capacity * sizeof(*history));
    if (!history) {
        return;
    }

    uint32_t path_metrics[NUM_STATES];
    uint32_t old_metrics[NUM_STATES];
    
    // Initiera felkostnader (tillstånd 0 är garanterad start)
    for (int i = 1; i < NUM_STATES; i++) {
        path_metrics[i] = 1e6;
    }
    path_metrics[0] = 0;

    int ch;

    // 1. Framåtpass (Forward Pass / Add-Compare-Select)
    while ((ch = fgetc(in)) != EOF) {
        for (int bit_pair = 3; bit_pair >= 0; bit_pair--) {
            uint8_t b0 = (ch >> (bit_pair * 2 + 1)) & 1;
            uint8_t b1 = (ch >> (bit_pair * 2)) & 1;

            // Expandera historikbufferten dynamiskt om det behövs
            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;
            }

            memcpy(old_metrics, path_metrics, sizeof(path_metrics));

            // Gå igenom alla möjliga måltillstånd
            for (int state = 0; state < NUM_STATES; state++) {
                int prev0 = state >> 1;
                int prev1 = prev0 | (NUM_STATES >> 1);
                uint8_t input_bit = state & 1;

                // Väg 0
                uint8_t word0 = (prev0 << 1) | input_bit;
                uint8_t exp0_0 = parity(word0 & POLY_G0);
                uint8_t exp0_1 = parity(word0 & POLY_G1);
                uint32_t cost0 = (exp0_0 != b0) + (exp0_1 != b1);

                // Väg 1
                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);

                uint32_t m0 = old_metrics[prev0] + cost0;
                uint32_t m1 = old_metrics[prev1] + cost1;

                // Välj den bästa överlevande vägen
                if (m0 <= m1) {
                    path_metrics[state] = m0;
                    history[step][state] = 0;
                } else {
                    path_metrics[state] = m1;
                    history[step][state] = 1;
                }
            }
            step++;
        }
    }

    if (step == 0) {
        free(history);
        return;
    }

    // 2. Hitta bästa sluttillstånd (lägsta 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;
        }
    }

    // 3. Bakåtspårning (Traceback)
    uint8_t *decoded_bits = malloc(step);
    if (!decoded_bits) {
        free(history);
        return;
    }

    for (int s = (int)step - 1; s >= 0; s--) {
        decoded_bits[s] = curr_state & 1;
        int decision = history[s][curr_state];
        int prev0 = curr_state >> 1;
        
        if (decision == 0) {
            curr_state = prev0;
        } else {
            curr_state = prev0 | (NUM_STATES >> 1);
        }
    }

    // 4. Återpacka till bytes och skriv ut till strömmen
    uint8_t current_byte = 0;
    for (size_t i = 0; i < step; i++) {
        current_byte = (current_byte << 1) | decoded_bits[i];
        if ((i + 1) % 8 == 0) {
            fputc(current_byte, out);
            current_byte = 0;
        }
    }

    free(decoded_bits);
    free(history);
}

// ============================================================================
// INTERNT TEST (-t)
// ============================================================================
void run_test(void) {
    printf("[*] Startar internt test (64 bytes randomiserat data)...\n");
    
    uint8_t original_data[64];
    srand((unsigned int)time(NULL));
    for (int i = 0; i < 64; i++) {
        original_data[i] = rand() % 256;
    }

    FILE *orig_stream = tmpfile();
    FILE *encoded_stream = tmpfile();
    FILE *decoded_stream = tmpfile();

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

    fwrite(original_data, 1, 64, orig_stream);
    rewind(orig_stream);

    encode(orig_stream, encoded_stream);
    rewind(encoded_stream);

    decode(encoded_stream, decoded_stream);
    rewind(decoded_stream);

    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 {
        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;
            }
        }
    }

    fclose(orig_stream);
    fclose(encoded_stream);
    fclose(decoded_stream);

    if (success) {
        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");
}

int main(int argc, char *argv[]) {
    if (argc != 2) {
        print_help(argv[0]);
        return EXIT_FAILURE;
    }

    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();
    } else {
        fprintf(stderr, "Fel: Okänd parameter '%s'\n", argv[1]);
        print_help(argv[0]);
        return EXIT_FAILURE;
    }

    return EXIT_SUCCESS;
}