Codici auto-correttori

Giulio Cesare Barozzi
CIRAM - Università di Bologna

[Revisione:luglio 2001. http://eulero.ing.unibo.it/~barozzi/MI2]

Richard W. Hamming, 1915-1998

Clear["Global`"]

Introduciamo la matrice di Hamming:

H = {{0,0,0,1,1,1,1},
	 {0,1,1,0,0,1,1},
	 {1,0,1,0,1,0,1}}; 
MatrixForm[%]

Le colonne della matrice di Hamming rappresentano (in base 2) i numeri da 1 a 7; dunque sono tutte le terne non nulle di ((GF(2))^3). Ricordiamo che la somma di due vettori non nulli di uno spazio vettoriale ((GF(2))^n)è ancora non nulla, a meno che tali vettori non siano tra loro uguali: x + y = 0  ⇒ x = y. Dunque non soltanto le colonne della matrice H sono non nulle, ma anche i vettori che si ottengono sommando due colonne sono diversi dal vettore nullo. 
Consideriamo due vettori di ((GF(2))^7), che indichiamo  x  e  y :

x = {1,0,1,0,1,0,1}; y = {1,0,0,1,1,1,1};

Verifichiamo che x  appartiene al nucleo ( = spazio nullo = kernel ) della matrice H , mentre y non gli appartiene:

{Mod[H.x,2], Mod[H.y, 2]}

I vettori che appartengono al nucleo della matrice di Hamming (si tratta di un sottospazio vettoriale di ((GF(2))^7), costituito da 16 elementi) sono le parole del codice, dunque  x  è una parola-codice, y non lo è. 
Vogliamo studiare più da vicino il nucleo della matrice H. Selezioniamo tre colonne di tale matrice che formino un minore di ordine 3 diverso da 0; ad esempio le ultime tre colonne. Procuriamoci una funzione che preleva i primi quattro elementi di una lista:

primi4[lista_]:= Take[lista, 4]

H1 = Map[primi4, H]; 
MatrixForm[%]

Una funzione che preleva gli ultimi tre elementi di una lista:

ultimi3[lista_] := Take[lista, -3]

H2 = Map[ultimi3, H];
MatrixForm[%]

Verifichiamo che H2 è invertibile:

Det[%]

Calcoliamo l'inversa di H2 :

Inverse[H2]//MatrixForm

... e riduciamola modulo 2:

invH2 = Mod[%, 2];
MatrixForm[%]

Moltiplichiamo l'inversa di  H2  per  H1 :

matKer = Mod[invH2.H1, 2];
MatrixForm[%]

Per ogni vettore x , ricaviamo  x5, x6, x7  in funzione di  x1, x2, x3, x4:

{x5, x6, x7} = matKer.{x1, x2, x3, x4}; 
MatrixForm[%]

Una generica parola del codice si scrive

c = Join[{x1, x2, x3, x4}, {x5, x6, x7}]; 
MatrixForm[%]

Elenchiamo tutte le parole del codice:

Table[c,{x1, 0,1},{x2, 0,1},{x3, 0,1},{x4, 0,1}];

codice = Mod[Flatten[%, 3], 2];
MatrixForm[%]

Riprendiamo in considerazione un vettore non appartenente a nucleo di H, ad esempio

x ={1,0,0,1,1,1,1};

Verifica:

Mod[H.x, 2]

... oppure

MemberQ[x, codice]

Ma allora il vettore a tre componenti H.x  coincide necessariamente con una delle colonne di H : da un esame diretto si riconosce che si tratta della prima colonna.
In generale, se  H.x  non è  il vettore nullo, come fare per trovare l'indice della colonna di  H  che fornisce tale vettore? Basta interpretare le tre componenti di H.x  come le cifre in base 2 di tale indice. Introduciamo una terna con le potenze di 2:

bin = {4, 2, 1};

Basta allora calcolare  Mod[H.x, 2].bin  per avere l'indice richiesto. Ad esempio, se x è il vettore precedentemente considerato:

Mod[H.x, 2].bin

In effetti

Mod[H.x,2]

è la prima colonna della matrice H.

Procuriamoci i vettori della base canonica di ((GF(2))^7). Cominciamo dalla matrice unità di dimensione 7:

e = IdentityMatrix[7]; e//MatrixForm

Le righe di tale matrice ci forniscono i vettori in questione. Consideriamo ancora il vettore  x  di prima. Visto che H.x coincide con la prima colonna di  H e questa non è altro che  H.e[[1]], andiamo a sommare (modulo 2) x  con e[[1]], cioè invertiamo il primo bit di x :

Mod[x + e[[1]], 2]

Dico che il risultato è una parola del codice:

Mod[H.%, 2]

Costruiamo un criterio per la correzione dei vettori in ingresso. Se la parola in ingresso (è la parola ricevuta) appartiene al codice, essa viene restituita in uscita, altrimenti viene restituita in uscita la parola (appartenente al codice) che si ottiene "correggendo" (cioè invertendo un bit della parola in ingresso). Precisamente il bit corrispondente alla posizione della colonna di H che coincide con  H.x

corr[x_] :=
	Module[{bin = {4, 2, 1}, e = IdentityMatrix[7], p},
		p = Mod[H.x, 2].bin;
		If[p == 0, x, Mod[x + e[[p]], 2] ]
	]

Proviamo con il vettore  x di prima:

{x, corr[x]}

Introduciamo in  ((GF(2))^7)la distanza di Hamming (si tratta di una distanza di indice 1): essa conta in quante posizioni i bit del vettore x differiscono dai corrispondenti bit del vettore y :

distH[x_, y_] := Apply[Plus, Mod[x + y, 2]]

Calcoliamo le mutue distanze tra gli elementi del codice:

Table[
        distH[codice[[i]],  
          codice[[j]]], {i, 16}, {j, 
          16}] // MatrixForm

Come si vede, la minima distanza tra elementi distinti del codice è 3 (in gergo DMIN = 3). Il risultato non è casuale. Se  x   è un vettore del codice , i vettori a distanza 1 da esso si scrivono  
	x + e[[h]] , 
quelli a distanza 2 si scrivono 
	x + e[[h]] + e[[k]] , 
dove 1 ≤ h ≠ k ≤ 7. Dunque l'applicazione della matrice  H a tali vettori fornisce la colonna di indice h (oppure la somma tra le colonne di indice h e k) della stessa matrice; in ogni caso un vettore non nullo, e dunque i vettori considerati non appartengono al codice.
Per ogni x  non appartenente al codice la funzione corr[x] fornisce l'unico elemento del codice che dista 1 da esso (l'elemento corretto è stato ottenuto dall'elemento in ingresso correggendo un solo bit). Tale elemento del codice è unico: se ce ne fossero due, essi disterebbero tra loro al più 2 (in virtù della diseguaglianza triangolare), in contraddizione con quanto abbiamo appena mostrato a proposito delle mutue distanze tra elementi del codice. 
In termini più formali: la funzione corr ha come dominio ((GF(2))^7) e come immagine codice, sottospazio di \(TraditionalForm\`dimensione\ 4\ di \ \(GF(2)\)\^7\).

Data una parola del codice, possiamo elencare tutti i vettori di ((GF(2))^7)a distanza 1 da esso. Costruiamo dapprima tutto ((GF(2))^7):

GF7  =  
        Flatten[Table[{x1, x2, x3, x4, x5, x6, x7}, 
{x1, 
              0, 1}, {x2, 0, 1}, {x3, 0, 1}, {x4, 0, 1}, {x5, 0, 1}, {x6, 0, 
              1}, {x7, 0, 1}], 6];

e poi selezioniamo gli elementi a distanza 1 da un vettore del codice, ad esempio:

Select[
      GF7,  (distH[
              codice[[10]], #] == 
            1) &]

Se viene ricevuto un elemento del codice, o esso coincide con l'elemento trasmesso (cioè non si sono avuti errori in trasmissione) oppure si sono avuti almeno 3 errori. Ammettendo che la probabilità di inversione di un bit in trasmissione sia p, e supponendo che gli eventuali errori siano eventi tra loro indipendenti, mettiamo a confronto la probabilità di 0 errori con la probabilità di 3 errori:

prob0 = (1 - p)^7; prob3 = Binomial[7, 3] (1-p)^4 p^3;

rapp = prob0/prob3

Esaminiamo il rapporto ottenuto in funzione di p :

Plot[rapp, {p, .001, .5}, 
PlotRange -> {0, 3},Frame -> True, GridLines -> {Automatic,{1., 2.}}];

Come si vede tale rapporto è nettamente maggiore di 1 per tutti i valori ragionevolmente piccoli di p. Ad esempio:

rapp/.p -> 0.000001

Calcoliamo il valore di p in corrispondenza del quale il rapporto vale 1:

rapp - 1;

Together[%]

Solve[Numerator[%] == 0]

Il primo valore fornito è reale:

First[%//N]

Se durante la trasmissione è avvenuto un errore, l'elemento ricevuto viene effetivamente corretto (cioè l'elemento restituito dalla funzione corr è quello effetivamente trasmesso), se invece avvengono due errori l'elemento restituito dalla funzione corr è un elemento del codice  diverso da quello trasmesso. Il codice che stiamo esaminando è dunque un SECC, cioè Single Error Correcting Code.
Calcoliamo la probabilità di 1 errore e di 2 errori:

prob1 = 7 (1-p)^6 p;

prob2 = 21 (1-p)^5 p^2;

Confrontiamo prob1 con prob2:

rapp = prob1/prob2;

Plot[rapp, {p,.01, 0.5}, 
PlotRange -> {0, 2}, Frame -> True, GridLines -> {Automatic,{1., 2.}}];

Together[rapp-1]

Il valore che annulla il numeratore è chiaramente 0.25 = 1/4. Infine confrontiamo la probabilità di 1 errore e la probabilità di 3 errori:

rapp = prob1/prob3;

Plot[rapp, {p, .01, .5},
PlotRange -> {0,2}, Frame -> True, GridLines -> {Automatic,{1., 2.}}];

Together[rapp-1]

Solve[Numerator[%] == 0]

Last[%]//N

Visualizzazione

Costruiamo una funzione che trasforma ogni vettore di ((GF(2))^(7 ))in una poligonale con vertici a coordinate intere, contenuta nel primo quadrante, con il seguente criterio: leggendo i bit da sinistra a destra, ogni bit uguale e 0 determina uno spostamento di un passo a destra, ogni bit uguale a uno determina uno spostamento verso l'alto. Ad esempio sia

x  =  {1, 0, 0, 1, 1, 0, 1};

Il vettore che ha per componenti i complementi a 1 delle componenti di x:

y = Mod[x + 1, 2]

Calcoliamo le "somme parziali" a partire da x e da y

ord = FoldList[Plus, 0, x]

asc  =  FoldList[Plus, 0, y]

... poi i punti corrispondenti

punti = Transpose[{asc, ord}]

Calcoliamo la poligonale:

traccia[x_] := 
      Module[{asc,  ord, y}, 
ord =  
          FoldList[Plus, 0, x]; 
y = Mod[x + 1, 2]; 
        asc  =  FoldList[Plus, 0, y]; 
punti = 
          Transpose[{asc, ord}]; 
punti]

Ad esempio:

traccia[
      codice[[11]] ]

Tracciamo la poligonale:

grafPoli[x_]  :=  
        ListPlot[traccia[x], PlotJoined -> True, 
          AspectRatio -> 
            Automatic, 
PlotRange -> {{0, 7}, {0, 7}}, 
          PlotStyle -> {RGBColor[1, 0, 0], 
              Thickness[0.01]}, 
DisplayFunction -> 
            Identity];

Tracciamo un reticolo di riferimento:

retUno  =  Table[{{0, k}, {7 - k, k}}, {k, 1, 7}];

retDue  =  Map[RotateLeft, %, {2}];

reticolo  = Join[retUno, retDue];

grafReticolo  = Graphics[Map[Line,  reticolo]];

Tutto insieme, per visualizzare il 13-esimo elemento del codice :

Show[{grafReticolo, 
          grafPoli[ codice[[13]] ] }, 
        DisplayFunction -> $DisplayFunction, 
        AspectRatio -> Automatic, Axes -> True, 
        Ticks -> {Range[7], Range[7]}];

Per rappresentare le 16 parole del codice basterebbe attivare il comando seguente (meglio non farlo, se si dispone di poca memoria):

Mettiamo a confronto i grafici di un vettore non appartenente al codice e del vettore che rappresenta la sua "correzione":

{x,  corr[x]}

Show[{grafReticolo, grafPoli[ x ] }, 
        DisplayFunction -> Identity, AspectRatio -> Automatic, 
        Axes -> True, Ticks -> {Range[7], Range[7]}];

Show[{grafReticolo, grafPoli[corr[x ]] }, 
        DisplayFunction -> Identity, AspectRatio -> Automatic, 
        Axes -> True, Ticks -> {Range[7], Range[7]}];

Show[GraphicsArray[{%%, %}], 
        DisplayFunction -> $DisplayFunction];

Come si vede la"correzione" ha modificato il bit numero 7; in generale l'indice del bit da modificare è dato da

p  = Mod[H . x, 2] . bin

come è mostrato nel corpo della funzione corr. Cerchiamo di evidenziare con un diverso colore il lato della poligonale che corrisponde alla correzione. Il segmento da evidenziare è:

Take[traccia[corr[x]], {p, p + 1}]

Definiamo una funzione che individui il segmento in cui la poligonale che rappresenta corr[x] si distacca da quella che rappresenta x :

segmento[x_, p_] :=Line[Take[traccia[x], {p, p + 1}]];

Tutto insieme:

graf[
          x_] := 
Module[{y =  corr[x], 
            p  = Mod[H . x, 2] . bin }, 
  If[
            y == 
              x, 
Show[{grafReticolo, grafPoli[ x ] }, 
              DisplayFunction -> $DisplayFunction, 
              AspectRatio -> Automatic, Axes -> True, 
              Ticks -> {Range[7], 
                  Range[7]}], 
Show[{grafReticolo, 
                grafPoli[x], grafPoli[y], 
                Graphics[{RGBColor[0, 0, 1], Thickness[0.01], 
                    segmento[y, p]}] }, 
              DisplayFunction -> $DisplayFunction, 
              AspectRatio -> Automatic, Axes -> True, 
              Ticks -> {Range[7], 
                  Range[7]}]
]
];

Ad esempio:

graf[x];

Sperimentazione:

graf[Table[Random[Integer, {0, 1}], {7}]  ];

La poligonale contenente il segmento blu (se esiste) rappresenta corr[x], altrimenti viene mostrata soltanto la poligonale corrispondente a x.