Announcement

Collapse

Forum Guidelines

This forum is for finished source code that is working properly. If you have questions about this or any other source code, please post it in one of the Discussion Forums, not here.
See more
See less

Reed Solomon Error Checking

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Reed Solomon Error Checking

    Hi everybody,

    The following is an adaptation of C code that handles Reed Solomon error correction for data being transmitted. I don't really understand what is going on under the hood (duh) but did the translation from C to PowerBasic. I thought it would be a cool exercise and something fun to do.

    where is this used for ? Well on usenet you may see lots of PAR2 file which are parity archives and QuickPar uses the algorithm internally to correct wrong packets. You could use it for UDP protocols f.i.

    Cheers

    Code:
    ' ======================================================================================================
    ' Reed Solomon Codes
    ' ======================================================================================================
    '/*             rs.c        */
    '/* This program is an encoder/decoder for Reed-Solomon codes. Encoding is in
    '   systematic form, decoding via the Berlekamp iterative algorithm.
    '   In the present form , the constants mm, nn, tt, and kk=nn-2tt must be
    '   specified  (the double letters are used simply to avoid clashes with
    '   other n,k,t used in other programs into which this was incorporated!)
    '   Also, the irreducible polynomial used to generate GF(2**mm) must also be
    '   entered -- these can be found in Lin and Costello, and also Clark and Cain.
    
    '   The representation of the elements of GF(2**m) is either in index form,
    '   where the number is the power of the primitive element alpha, which is
    '   convenient for multiplication (add the powers modulo 2**m-1) or in
    '   polynomial form, where the bits represent the coefficients of the
    '   polynomial representation of the number, which is the most convenient form
    '   for addition.  The two forms are swapped between via lookup tables.
    '   This leads to fairly messy looking expressions, but unfortunately, there
    '   is no easy alternative when working with Galois arithmetic.
    
    '   The code is not written in the most elegant way, but to the best
    '   of my knowledge, (no absolute guarantees!), it works.
    '   However, when including it into a simulation program, you may want to do
    '   some conversion of global variables (used here because I am lazy!) to
    '   local variables where appropriate, and passing parameters (eg array
    '   addresses) to the functions  may be a sensible move to reduce the number
    '   of global variables and thus decrease the chance of a bug being introduced.
    
    '   This program does not handle erasures at present, but should not be hard
    '   to adapt to do this, as it is just an adjustment to the Berlekamp-Massey
    '   algorithm. It also does not attempt to decode past the BCH bound -- see
    '   Blahut "Theory and practice of error control codes" for how to do this.
    
    '              Simon Rockliff, University of Adelaide   21/9/89
    
    '   26/6/91 Slight modifications to remove a compiler dependent bug which hadn't
    '           previously surfaced. A few extra comments added for clarity.
    '           Appears to all work fine, ready for posting to net!
    
    '                  Notice
    '                 --------
    '   This program may be freely modified and/or given to whoever wants it.
    '   A condition of such distribution is that the author's contribution be
    '   acknowledged by his name being left in the comments heading the program,
    '   however no responsibility is accepted for any financial or other loss which
    '   may result from some unforseen errors or malfunctioning of the program
    '   during use.
    '                                 Simon Rockliff, 26th June 1991
    '
    '   ADAPTED TO POWERBASIC (TM) BY STEVEN PRINGELS
    '*/
    
    
    %mm = 4
    %nn = 15
    %tt = 3
    %kk = 9
    ' ======================================================================================================
    Function initReedSolomon() As Long
      Dim pp      (%mm + 1)       As Global Integer
      Dim alpha_to(%nn + 1)       As Global Integer
      Dim index_of(%nn + 1)       As Global Integer
      Dim gg      (%nn - %kk + 1) As Global integer
      Dim recd    (%nn)           As Global integer
      Dim data_int(%kk)           As Global Integer
      Dim bb      (%nn - %kk)     As Global integer
    
      pp(0) = 1
      pp(1) = 1
      pp(2) = 0
      pp(3) = 0
      pp(4) = 1
    
    End Function
    
    Function generate_gf() As Long
      Local Mask    As Integer
      Local i       As Integer
      Local Temp    As integer
    
      mask          = 1
      alpha_to(%mm) = 0
    
      For i = 0 To %mm - 1
    
        alpha_to(i) = mask
    
         index_of(alpha_to(i)) = i
         If pp(i) <> 0 Then
           alpha_to(%mm) = alpha_to(%mm) Xor mask
         End if
         Shift Left mask, 1
    
      Next
    
      index_of(alpha_to(%mm)) = %mm
      Shift Right mask, 1
      for i = %mm + 1 To %nn - 1
         if alpha_to(i-1) >= mask Then
    
            Temp = (alpha_to(i-1) Xor mask)
    
            Shift Left Temp, 1
    
            alpha_to(i) = alpha_to(%mm) Xor Temp
    
         else
            Temp = alpha_to(i-1)
            Shift Left Temp, 1
            alpha_to(i) = Temp
         End iF
         index_of(alpha_to(i)) = i
      Next
      index_of(0) = -1
    
    End Function
    
    Function gen_poly() As Long
    
    '/* Obtain the generator polynomial of the tt-error correcting, length
    '  nn=(2**mm -1) Reed Solomon code  from the product of (X+alpha**i), i=1..2*tt
    '*/
    
       Local i, j As integer
    
       gg(0) = 2    ' /* primitive element alpha = 2  for GF(2**mm)  */
       gg(1) = 1    ' /* g(x) = (X+alpha) initially */
    
       for i= 2 To (%nn - %kk)
          gg(i) = 1
          for j = i - 1 To 1 Step -1
              if gg(j) <> 0 Then
                 gg(j) = gg(j-1) Xor alpha_to((index_of(gg(j)) + i ) mod %nn)
              else
                 gg(j) = gg(j-1)
              End If
          next
          gg(0) = alpha_to((index_of(gg(0)) + i ) mod %nn) '     /* gg[0] can never be zero */
       Next
       ' /* convert gg[] to index form for quicker encoding */
       For i=0 To (%nn-%kk)
           gg(i) = index_of(gg(i))
       Next
    
    End Function
    
    Function encode_rs() As Long
    
    '/* take the string of symbols in data[i], i=0..(k-1) and encode systematically
    '   to produce 2*tt parity symbols in bb[0]..bb[2*tt-1]
    '   data[] is input and bb[] is output in polynomial form.
    '   Encoding is done by using a feedback shift register with appropriate
    '   connections specified by the elements of gg[], which was generated above.
    '   Codeword is   c(X) = data(X)*X**(nn-kk)+ b(X)          */
    
       Local i, j     As Integer
       Local feedback As integer
    
       for i=0 To (%nn - %kk) - 1
           bb(i) = 0
       next
    
       for i = %kk - 1 To 0 Step -1
           feedback = index_of(data_int(i) Xor bb(%nn - %kk - 1))
    
           if feedback <> -1 Then
              for j = %nn - %kk - 1 To 1 Step -1
                if (gg(j) <> -1) Then
                   bb(j) = bb(j-1) Xor alpha_to((gg(j) + feedback) mod %nn)
                else
                   bb(j) = bb(j-1)
                End If
              Next
              bb(0) = alpha_to((gg(0) + feedback) mod %nn)
           else
              for j = %nn - %kk - 1 To 1 Step -1
                  bb(j) = bb(j-1)
              Next
              bb(0) = 0
           End If
       Next
    
    End Function
    
    Function decode_rs() As Long
    
    '/* assume we have received bits grouped into mm-bit symbols in recd[i],
    '   i=0..(nn-1),  and recd[i] is index form (ie as powers of alpha).
    '   We first compute the 2*tt syndromes by substituting alpha**i into rec(X) and
    '   evaluating, storing the syndromes in s[i], i=1..2tt (leave s[0] zero) .
    '   Then we use the Berlekamp iteration to find the error location polynomial
    '   elp[i].   If the degree of the elp is >tt, we cannot correct all the errors
    '   and hence just put out the information symbols uncorrected. If the degree of
    '   elp is <=tt, we substitute alpha**i , i=1..n into the elp to get the roots,
    '   hence the inverse roots, the error location numbers. If the number of errors
    '   located does not equal the degree of the elp, we have more than tt errors
    '   and cannot correct them.  Otherwise, we then solve for the error value at
    '   the error location and correct the error.  The procedure is that found in
    '   Lin and Costello. For the cases where the number of errors is known to be too
    '   large to correct, the information symbols as received are output (the
    '   advantage of systematic encoding is that hopefully some of the information
    '   symbols will be okay and that if we are in luck, the errors are in the
    '   parity part of the transmitted codeword).  Of course, these insoluble cases
    '   can be returned as error flags to the calling routine if desired.   */
    
       Local i,j,u,q As Integer
    
       Dim   elp (%nn - %kk + 2, %nn - %kk) As Integer
       Dim   d   (%nn - %kk + 2)            As integer
       Dim   l   (%nn - %kk + 2)            As Integer
       Dim   u_lu(%nn - %kk + 2)            As Integer
       Dim   s   (%nn - %kk + 1)            As Integer
    
       'int elp[nn-kk+2][nn-kk], d[nn-kk+2], l[nn-kk+2], u_lu[nn-kk+2], s[nn-kk+1] ;
       Local Count     As Integer
       Local syn_error As integer
    
       Dim root(%tt)    As Integer
       Dim loc_int(%tt) As Integer
       Dim z(%tt + 1)   As integer
       Dim err_int(%nn) As Integer
       Dim reg(%tt + 1) As Integer
    
    
      ' /* first form the syndromes */
       For i = 1 To %nn - %kk
           s(i) = 0
    
           for j = 0 To %nn - 1
               if recd(j) <> -1 Then
                  s(i) = s(i) Xor alpha_to((recd(j) + i * j ) mod %nn) '    /* recd[j] in index form */
               End If
           next
    
           ' /* convert syndrome from polynomial form to index form  */
           If (s(i) <> 0) Then
              syn_error = 1
           End If
           s(i) = index_of(s(i))
    
       Next
    
       if syn_error Then ' /* if errors, try and correct */
    
    '/* compute the error location polynomial via the Berlekamp iterative algorithm,
    '   following the terminology of Lin and Costello :   d[u] is the 'mu'th
    '   discrepancy, where u='mu'+1 and 'mu' (the Greek letter!) is the step number
    '   ranging from -1 to 2*tt (see L&C),  l[u] is the
    '   degree of the elp at that step, and u_l[u] is the difference between the
    '   step number and the degree of the elp.
    '*/
    '/* initialise table entries */
    
          d(0)     = 0          ' /* index form */
          d(1)     = s(1)       ' /* index form */
          elp(0,0) = 0          ' /* index form */
          elp(1,0) = 1          ' /* polynomial form */
    
          for i=1 To (%nn - %kk) - 1
             elp(0, i) = -1 '  /* index form */
             elp(1, i) = 0  '  /* polynomial form */
          Next
    
          l(0) = 0
          l(1) = 0
          u_lu(0) = -1
          u_lu(1) = 0
          u = 0
    
          do
            Incr u
            if d(u) = -1 Then
               l(u+1) = l(u)
               for i = 0 To l(u)
                   elp(u+1 , i) = elp(u, i)
                   elp(u, i) = index_of(elp(u,i))
               Next
            else
            ' /* search for words with greatest u_lu[q] for which d[q]!=0 */
               q = u - 1
               while ((d(q) = -1) And (q > 0 ))
                 Decr q
               Wend
    
               ' /* have found first non-zero d[q]  */
               if (q > 0) Then
                   j=q
                   do
                     Decr j
                     if ((d(j) <> -1 ) And (u_lu(q) < u_lu(j))) Then
                         q = j
                     End If
                   Loop while j > 0
               End If
    
               ' /* have now found q such that d[u]!=0 and u_lu[q] is maximum */
               ' /* store degree of new elp polynomial */
               if ( l(u) > l(q) + u - q ) Then
                  l(u+1) = l(u)
               else
                  l(u+1) = l(q) + u - q
               End If
    
               ' /* form new elp(x) */
               for i = 0 To (%nn - %kk) - 1
                   elp(u+1, i) = 0
               Next
    
               for i = 0 To l(q)
                   if (elp(q, i) <> -1) Then
                      elp(u + 1, i + u - q) = alpha_to((d(u) + %nn - d(q) + elp(q, i)) Mod %nn)
                   End If
               Next
    
               For i = 0 To l(u)
                   elp(u+1, i) = elp(u + 1, i) Xor elp(u, i)
                   elp(u, i) = index_of(elp(u, i)) '  /*convert old elp value to index*/
               Next
    
            End If
    
            u_lu( u + 1) = u - l(u+1)
    
            Rem ./* form (u+1)th discrepancy */
            if (u < %nn - %kk) Then    '/* no discrepancy computed on last iteration */
               if (s(u+1) <> -1) Then
                   d(u+1) = alpha_to(s(u+1))
               else
                   d(u+1) = 0
               End If
               for i = 1 To l(u + 1)
                   if ((s( u + 1 - i) <> -1) And (elp(u+1, i) <> 0 ) ) Then
                      d( u + 1 ) = d(u + 1) Xor alpha_to((s(u + 1 - i) + index_of(elp( u + 1, i))) mod %nn)
                   End If
               Next
               d( u + 1 ) = index_of(d(u+1)) ' /* put d[u+1] into index form */
            End If
    
          Loop while (( u < %nn - %kk) And (l(u+1) <= %tt))
    
          Incr u
    
          if (l(u) <= %tt) Then       '  /* can correct error */
             '/* put elp into index form */
    
             for i=0 To l(u)
                 elp(u, i) = index_of(elp(u, i))
             Next
             '/* find roots of the error location polynomial */
    
             for i=1 To i <= l(u)
                 reg(i) = elp(u, i)
             Next
    
             count = 0
    
             For i=1 To %nn
                 q = 1
                 for j=1 To l(u)
                  if (reg(j) <> -1) Then
                     reg(j) = (reg(j)+j) mod %nn
                     q = q Xor alpha_to(reg(j))
                  End If
                 Next
                 if Not q Then        ' /* store root and error location number indices */
                    root(count) = i
                    loc_int(count) = %nn - i
                    Incr count
                 End If
             Next
             if (count = l(u)) Then   ' /* no. roots = degree of elp hence <= tt errors */
    
                ' /* form polynomial z(x) */
                for i=1 To l(u)       ' /* Z[0] = 1 always - do not need */
                    If ((s(i) <> -1) And (elp(u, i) <> -1)) Then
                       z(i) = alpha_to(s(i)) Xor alpha_to(elp(u, i))
                    Else
                       if ((s(i) <> -1) And (elp(u, i) = -1)) Then
                          z(i) = alpha_to(s(i))
                       else
                           if ((s(i) = -1) And (elp(u, i) <> -1)) Then
                              z(i) = alpha_to(elp(u, i))
                           else
                              z(i) = 0
                           End If
                       End If
                    End If
                    For j = 1 To i - 1
                       if ((s(j) <> -1) And (elp(u, i-j) <> -1)) Then
                          z(i) = z(i) Xor alpha_to((elp(u, i-j) + s(j)) mod %nn)
                       End If
                    next
                    z(i) = index_of(z(i)) '/* put into index form */
                Next
    
               ' /* evaluate errors at locations given by error location numbers loc[i] */
               for i=0 To %nn - 1
                   err_int(i) = 0
                   if (recd(i) <> -1) Then        ' /* convert recd[] to polynomial form */
                       recd(i) = alpha_to(recd(i))
                   else
                       recd(i) = 0
                   End If
               Next
    
               for i=0 To l(u) - 1        ' /* compute numerator of error term first */
                  err_int(loc_int(i)) = 1 ' /* accounts for z[0] */
                  for j=1 To l(u)
                      if (z(j) <> -1) Then
                         err_int(loc_int(i)) = err_int(loc_int(i)) Xor alpha_to((z(j) + j * root(i)) mod %nn)
                      End If
                  Next
                  if (err_int(loc_int(i)) <> 0 ) Then
                      err_int(loc_int(i)) = index_of(err_int(loc_int(i)))
                      q = 0 ' /* form denominator of error term */
                      for j=0 To l(u) - 1
                         if (j <> i) Then
                            q = q + index_of(1 xor alpha_to((loc_int(j) + root(i)) mod %nn))
                         End If
                      next
                      q = q mod %nn
                      err_int(loc_int(i)) = alpha_to((err_int(loc_int(i)) - q + %nn) mod %nn)
    
                      recd(loc_int(i)) = recd(loc_int(i)) xor err_int(loc_int(i)) ' /*recd[i] must be in polynomial form */
                  End If
              Next
             Else                             ' /* no. roots != degree of elp => >tt errors and cannot solve */
               For i=0 To %nn - 1             ' /* could return error flag if desired */
                   if (recd(i) <> -1) Then    ' /* convert recd[] to polynomial form */
                      recd(i) = alpha_to(recd(i))
                   else
                      recd(i) = 0             ' /* just output received codeword as is */
                   End iF
               Next
             End if
         Else                                 ' /* elp has degree has degree >tt hence cannot solve */
           for i=0 To %nn - 1                 ' /* could return error flag if desired */
              if (recd(i) <> -1) Then         ' /* convert recd[] to polynomial form */
                 recd(i) = alpha_to(recd(i))
              else
                 recd(i) = 0                  ' /* just output received codeword as is */
              End If
           Next
         End If
       Else                                   ' /* no non-zero syndromes => no errors: output received codeword */
        For i=0 To %nn - 1
           if (recd(i) <> -1) Then            ' /* convert recd[] to polynomial form */
              recd(i) = alpha_to(recd(i))
           else
              recd(i) = 0
           End If
        next
      End If
    
    End Function
    
    Function pbmain()
    
      Call initReedSolomon()
    
      Local i As integer
    
      ' generate the Galois Field GF(2 ^ %mm)
    
      Call generate_gf()
    
      Stdout "Look-up tables for GF(2^" & Format$(%mm) & ")"
    
      for i = 0 To %nn
         StdOut Format$(i) & ", " & Format$(alpha_to(i)) & ", " & Format$(index_of(i))
      Next
      Print
    
      ' /* compute the generator polynomial for this RS code */
    
      Call gen_poly()
    
    
      ' for known data_int, stick a few numbers into a zero codeword. data_int is in
      ' polynomial form.
    
      for i=0 To %kk - 1
         data_int(i) = 0
      next
    
      ' /* for example, say we transmit the following message (nothing special!)
    
     data_int(0) = 8
     data_int(1) = 6
     data_int(2) = 8
     data_int(3) = 1
     data_int(4) = 2
     data_int(5) = 4
     data_int(6) = 15
     data_int(7) = 9
     data_int(8) = 9
    
    ' encode data_int[] to produce parity in bb[].  data_int input and parity output
    ' is in polynomial form
    '
     stdout "== Encoding =="
     Call encode_rs()
    
    '/* put the transmitted codeword, made up of data_int plus parity, in recd[] */
    
     for i=0 To (%nn - %kk) - 1
         recd(i) = bb(i)
    
         stdout "Encoding: " & Format$(i) & ", " & Format$(recd(i))
     next
    
     for i=0 To %kk - 1
         recd(i + %nn - %kk) = data_int(i)
         stdout "Recd =  : " & Format$(i) & ", " & Format$(recd(i))
     next
    
    '/* if you want to test the program, corrupt some of the elements of recd[]
    '   here. This can also be done easily in a debugger. */
    '/* Again, lets say that a middle element is changed */
    
     data_int(2) = 3
    
     for i = 0 To %nn - 1
         recd(i) = index_of(recd(i))         ' /* put recd[i] into index form */
     next
    
     ' /* decode recv[] */
     Call decode_rs()                        ' /* recd[] is returned in polynomial form */
    
     ' /* print out the relevant stuff - initial and decoded {parity and message} */
     Stdout "Results for Reed-Solomon code n, k, t = " & Format$(%nn) & "," & Format$(%kk) & "," & Format$(%tt)
     Stdout "i      data(i)       decoded"
     for i=0 To (%nn - %kk) - 1
         Stdout format$(i, "000") & "      " & Format$(bb(i)) & "            " & Format$(recd(i))
     next
     Stdout "======================================="
     for i = (%nn - %kk) To (%nn - 1)
         Stdout Format$(i, "000") & "      " & Format$(data_int(i - %nn + %kk)) & "             " & Format$(recd(i))
     next
    
    End Function
    So here we are, this is the end.
    But all that dies, is born again.
    - From The Ashes (In This Moment)
Working...
X