reedsolomondecoder.h

心电图小波零树压缩演算法的研究

/*
  *********************************************************************
  *                                                                   *
  *        Galois Field Arithmetic Library (version 0.0.1)            *
  *                                                                   *
  * Class: Reed-Solomon Decoder                                       *
  * Version: 0.0.1                                                    *
  * Author: Arash Partow - 2000                                       *
  * URL: http://www.partow.net/projects/galois/index.html             *
  *                                                                   *
  * Copyright Notice:                                                 *
  * Free use of this library is permitted under the guidelines and    *
  * in accordance with the most current version of the Common Public  *
  * License.                                                          *
  * http://www.opensource.org/licenses/cpl.php                        *
  *                                                                   *
  *********************************************************************
*/

#ifndef INCLUDE_REEDSOLOMONDECODER_H
#define INCLUDE_REEDSOLOMONDECODER_H

#include <string>
#include "GaloisField.h"
#include "GaloisFieldElement.h"
#include "GaloisFieldPolynomial.h"
#include "ReedSolomonBlock.h"


using namespace galois;

namespace reedsolomon
{

   class ReedSolomonDecoder
   {

      public:

       ReedSolomonDecoder(GaloisField* _gf, const unsigned int _gen_initial_index, const unsigned int& _code_length, const unsigned int& _fec_length)
       {
          gf                = _gf;
          gen_initial_index = _gen_initial_index;
          code_length       = _code_length;
          fec_length        = _fec_length;
          data_length       = code_length - fec_length;
          bit_length        = data_length * gf->pwr();

          GaloisFieldElement xgfe[2] = {
                                        GaloisFieldElement(gf, 0),
                                        GaloisFieldElement(gf, 1)
                                       };

          X = GaloisFieldPolynomial(gf,1,xgfe);
       };


      ~ReedSolomonDecoder(){}


       bool decode(ReedSolomonBlock& rsblock, const std::vector <unsigned int>& erasure)
       {
          errors_detected = 0;
          GaloisFieldPolynomial received(gf,code_length - 1);
          load_message(received,rsblock);

          GaloisFieldPolynomial syndrome;

          /*
             No errors have been found within given data block
             is syndrome result is 0
          */
          if (compute_syndrome(received,syndrome) == 0)
          {
             rsblock.errors_detected  = 0;
             rsblock.errors_corrected = 0;
             rsblock.unrecoverable    = false;
             return true;
          }

          std::vector <unsigned int> erasure_list = erasure;

          prepare_erasures(erasure_list);


          // find error + erasure locator polynomial
          GaloisFieldPolynomial LFSR  = BerlekampMassey(syndrome, erasure_list); // lambda

          // find roots of error locator polynomial
          std::vector < GaloisFieldElement > LFSR_roots;
          find_roots(LFSR, LFSR_roots);

          errors_detected = static_cast<unsigned int>(LFSR_roots.size());

          // error evaluator polynomial
          GaloisFieldPolynomial omega = compute_omega(LFSR,syndrome);

          GaloisFieldPolynomial LFSR_derivative = LFSR.derivative();

          for (unsigned int i = 0; i < LFSR_roots.size(); i++)
          {

             GaloisFieldElement current_root = LFSR_roots[i];
             GaloisFieldElement numerator    = omega(current_root) * (current_root ^ (gen_initial_index - 1));
             GaloisFieldElement denominator  = LFSR_derivative(current_root);

             if (denominator == 0)
             {
                rsblock.errors_detected = errors_detected;
                rsblock.unrecoverable = true;
                return false;
             }

             GFSymbol position = gf->index(LFSR_roots[i].inverse());

             received[position] -=  (numerator / denominator); // apply correction
             rsblock.data[code_length - position - 1] = static_cast<char>(received[position].poly());

             rsblock.errors_corrected++;
          }

          rsblock.errors_detected = errors_detected;
          return true;
       }


      private:

       void load_message(GaloisFieldPolynomial& received, const ReedSolomonBlock& rsblock)
       {
          /*
            Load message data into received polynomial
          */
          std::string message = "";
          if (!rsblock.single_seq)
            message = rsblock.data + rsblock.fec;
          else
            message = rsblock.data;

          for(unsigned int i = 0; i < code_length; i++)
          {
             received[code_length - i - 1] = static_cast<GFSymbol>(message[i]);
          }
       }


       void prepare_erasures(std::vector <unsigned int>& erasure)
       {
          /*
            Erasures positions are given as zero index positions, they must be reverted
            so that they correspond to the relative polynomial terms
          */
          for(unsigned int i = 0; i < erasure.size(); i++)
          {
             erasure[i] = code_length - erasure[i];
          }
       };


       int compute_syndrome(const GaloisFieldPolynomial& received, GaloisFieldPolynomial& syndrome)
       {
          int error_flag = 0;
          unsigned int syn_error_count = 0;

          // compute syndromes
          syndrome = GaloisFieldPolynomial(gf,fec_length - 1); // create syndrome polynomial of degree fec_length - 1

          for(unsigned int i = 0; i < fec_length; i++)
          {
             syndrome[i]  = received(GaloisFieldElement(gf,gf->alpha(gen_initial_index + i)));   // initial_index is verified.
             error_flag  |= syndrome[i].poly();

             if (syndrome[i].poly() != 0)
               syn_error_count++;
          }

          return error_flag;
       }


       void compute_gamma(GaloisFieldPolynomial& gamma, const std::vector <unsigned int>& erasure)
       {
          if (erasure.size() > 0)
          {
             gamma = (1 + (X * GaloisFieldElement(gf,erasure[0])));

             for (unsigned int e = 1; e < erasure.size(); e++)
             {
                gamma *= (1 + (X * GaloisFieldElement(gf,erasure[e])));
             }
          }
          else
            gamma = GaloisFieldPolynomial(gf,0);
       }


       GaloisFieldPolynomial compute_omega(const GaloisFieldPolynomial& lambda, const GaloisFieldPolynomial& syndrome/*, GaloisFieldPolynomial& gamma*/)
       {
          /*
            Compute error/erasure evaluator polynomial omega = (lambda * syndrome) mod x^fec_length
          */
          return ((lambda * syndrome) % fec_length);
       }


       void find_roots(const GaloisFieldPolynomial& poly, std::vector < GaloisFieldElement >& root_list)
       {
          /*
            Find roots of polynomial via Chien search method
          */
          root_list.clear();
          for(unsigned int i = 0; i <= code_length; i++)
          {
             if (poly(i) == 0)
             {
                root_list.push_back(GaloisFieldElement(gf,i));
                if (root_list.size() == poly.deg())
                  break;
             }
          }
       }


       void compute_discrepancy(GaloisFieldElement&    discrepancy,
                          const GaloisFieldPolynomial& lambda,
                          const GaloisFieldPolynomial& syndrome,
                                unsigned int                  r)
       {
          /*
            Compute the lambda discrepancy at the r'th loop of BMA
          */

          discrepancy = 0;
          for (unsigned int i = 0; i <= lambda.deg(); i++)
          {
             discrepancy += lambda[i] * syndrome[r - i];
          }
       }


       GaloisFieldPolynomial BerlekampMassey (GaloisFieldPolynomial&            syndrome,  // syndrome polynomial
                                              const std::vector <unsigned int>& erasure // number of elements in erasure locations
                                             )
       {

          /*
             Berlekamp-Massey Algorithm
             Identify the "minimal" length of the linear feed-back shift register (LFSR) which will
             generate the sequence equivalent to the syndrome.
          */
          GaloisFieldPolynomial lambda    = GaloisFieldElement(gf, 1); // LFSR Connection polynomial
          GaloisFieldPolynomial lambda_a  = lambda;
          GaloisFieldPolynomial lambda_b  = X;
          GaloisFieldElement discrepancy  = GaloisFieldElement(gf, 0);

          for(unsigned int round = 0; round < fec_length; round++)
          {
             compute_discrepancy(discrepancy,lambda,syndrome,round);

             if (discrepancy != 0)
             {
                // An new error has been detected
                lambda_a = lambda - discrepancy * lambda_b;

                if ((2 * lambda.deg()) < (round + 1))
                {
                   lambda_b = lambda / discrepancy;
                }

                lambda = lambda_a;

                lambda_b <<= 1;
             }
          }
          return lambda;
       }

      private:

       GaloisField*           gf;
       GaloisFieldPolynomial  X;
       unsigned int           gen_initial_index;
       unsigned int           code_length;
       unsigned int           fec_length;
       unsigned int           data_length;
       unsigned int           bit_length; // data length in bits

       //Temporary variable
       unsigned int           errors_detected;

   };

}

#endif