# MAPLE. April 1, 2020
# show how interleaving (unshifted) and Hamming code handle burst errors.
#
# MAPLE. March 16, 2020
#
# Simulate burst error channel with Markov chain.
with(LinearAlgebra):

 # We will corrupt N/7 codewords from the simplest Hamming code.
 N:=70000;
 interface(rtablesize=N+2); # Tell MAPLE we're willing to see this displayed
 pval := 0.4;               # probability of a bit error when in faulty mode
 T12val := 0.005;           # probability of moving from good mode to faulty mode at any point
 T21val := 0.1;             # probability of moving out of the faulty mode at any point
# Solve for steady state vector of this Markov chain.
assign( solve({  (1-T12val)*p1+T21val*p2 = p1, T12val*p1+(1-T21val)*p2=p2,p1+p2=1.0 }, {p1,p2} ) );


# This is the p for the binary entropy:
Hp := p2*pval;
# and here's the entropy
H := -Hp*log(Hp) - (1-Hp)*log(1.-Hp) ;


#############################################################
### Generating bursty signal
#############################################################
myrand := rand(0..10000);

# Generate "true" with probability p, else "false"
p := proc()   # Return true w prob pval;
 global pval;
 local x, ans;
 ans := false;
 x := myrand();
 if x > (1-pval)*10000 then ans := true; fi;
 RETURN(ans);
end;

# Generate "true" with probability T12val, else "false"
T12 := proc()   # Return true w prob T12val;
 global T12val;
 local x, ans;
 ans := false;
 x := myrand();
 if x > (1-T12val)*10000 then ans := true; fi;
 RETURN(ans);
end ;

# Generate "true" with probability T21val, else "false"
T21 := proc()   # Return true w prob T21val;
 global T21val;
 local x, ans;
 ans := false;
 x := myrand();
 if x > (1-T21val)*10000 then ans := true; fi;
 RETURN(ans);
end ;


# Initial count on number of bursts, initialize signal, first bit b = 0
nbursts := 0;
b := 0;
signal := Array(1..N):
printf(`%d`,b);
good := true:
for i to N 
 do 
    if i mod 100 = 0 then printf(`\n`); fi;   # Only 100 bits per line of output
    if not good and p() then b := 1; fi;
    # A transition into and out of faulty mode is signaled in the output by a blank space
    if     good and T12() then 
                          good := false; nbursts := nbursts+1;  
                          printf(` `); fi; # Now faulty, count this burst
    printf(`%d`,b);
    signal[i] := b; # Store it in this long string
    if not good and T21() then good := true;  printf(` `); fi; # Exit faulty mode.
    b := 0;
 od:

printf(`\n\n`);

#############################################################
#############################################################
  
printf(` Decoding Statistics: \n`);



#############################################################
##  Hamming code of length 7
#############################################################
# Encode with Hamming code H7.  But we will just use zero codewords.
Enc := proc(m) RETURN( [ 
m[1] + m[2]        + m[4] mod 2 , 
m[1]        + m[3] + m[4] mod 2 , 
m[1]                      mod 2 , 
       m[2] + m[3] + m[4] mod 2 , 
       m[2]               mod 2 , 
              m[3]        mod 2 , 
                     m[4] mod 2 ]);  end: # Not used in this file

# Syndrome calculation.
Synd := proc(r)  RETURN( [
                     r[4] + r[5] + r[6] + r[7] mod 2, 
       r[2] + r[3]               + r[6] + r[7] mod 2, 
r[1]        + r[3]        + r[5]        + r[7] mod 2 ]); end:


# Decode with Hamming code H7.  But we will just use zero codewords.
Dec :=  proc(r) 
local h, s , ans;
for h to 7 do ans[h]:=r[h]; od;
s := Synd(r);
h := 4*s[1]+2*s[2]+s[3]; # If nonzero assume this is the error location.
if h > 0 then ans[h] := 1 - ans[h]; fi;  # FLip one bit.
RETURN( [ seq(ans[h] , h=1..7 ) ] );
end:

# The original message should reside in positions 3,5,6,7 of the corrected received vector
Recover := proc( r ) local c; c := Dec( r ); RETURN( [ c[3],c[5],c[6],c[7] ] ); end:
#############################################################
#############################################################

# Hamming weight
wt := proc(z)  local h; RETURN( add( z[h], h=1..nops(z) ) ); end;


# Basic stats on the randomly generated signal
print(` Signal contains `,add(signal[h],h=1..N),` errors, grouped in `,nbursts,` bursts of noise. `);

# Initialize counts before decoding
nbad := 0:  nbiterrors := 0:
ncol := floor(N/7):  # Arrange signal as  7 x ncol array and decode each column
rmsg := Array(1..4*ncol): # Recovered message should be all zeros.
for i to ncol
 do
  x := Recover( [ seq(signal[i+h*ncol] ,h=0..6) ] );  # Pull out i-th column and recover
  if x <> [0,0,0,0] then nbad := nbad + 1; nbiterrors := nbiterrors + wt(x); fi;  # Record errors
  rmsg[4*i-3] := x[1]; 
  rmsg[4*i-2] := x[2]; 
  rmsg[4*i-1] := x[3]; 
  rmsg[4*i  ] := x[4]; 
 od:

# Basic stats on recovered message (length 4/7 of signal length)
print(` Result: `,nbad,` 4-tuples decoded incorrectly for a total of `,nbiterrors,` damaged bits `);

print(` We receive:`,rmsg);

quit