LDPC Performance

Jason
2013-04-05
2013-05-02
  • Jason
    Jason
    2013-04-05

    I used IT++ some time ago for simulating LDPC for DVB-S2. The performance for that standard was about the same as what was reported in the literature, but I tried the same for WiGig and found the performance wasn't very good. The waterfall region was about right, but the error floor was very high. Has anyone else tried using IT++ for short LDPC codes? I'll just paste the entire main.cpp below (please let me know if this is against the rules or etiquette of the site). The code to follow can be compiled with g++ - my full command was 'g++ -O3 -o sim -I /usr/local/include -L/usr/local/lib main.cpp -litpp'

    Note: I'm trying to move this from "Open Discussion" into "Help," but this is my first time posting on sourceforge, and I'm not sure how to do that.

    Jason

    #include <iostream>
    
    #include <vector>
    #include <algorithm>
    
    #include <math.h>
    
    #include <time.h>
    
    #include <itpp/itcomm.h>
    using namespace std;
    using namespace itpp;
    
    void ExpandWigig( LDPC_Code* result, LDPC_Parity* parityMatrix, unsigned int modeNum )
    {
      const unsigned int EXPANSION = 42;
    
      unsigned int mode_13_16[8][16] = {
        { 29, 30, 0, 8, 33, 22, 17, 4, 27, 28, 20, 27, 24, 23, -1, -1 },
        { 37, 31, 18, 23, 11, 21, 6, 20, 32, 9, 12, 29, 10, 0, 13, -1 },
        { 25, 22, 4, 34, 31, 3, 14, 15, 4, 2, 14, 18, 13, 13, 22, 24 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }
      };
      unsigned int mode_3_4[8][16] = {
        { 35, 19, 41, 22, 40, 41, 39, 6, 28, 18, 17, 3, 28, -1, -1, -1 },
        { 29, 30, 0, 8, 33, 22, 17, 4, 27, 28, 20, 27, 24, 23, -1, -1 },
        { 37, 31, 18, 23, 11, 21, 6, 20, 32, 9, 12, 29, -1, 0, 13, -1 },
        { 25, 22, 4, 34, 31, 3, 14, 15, 4, -1, 14, 18, 13, 13, 22, 24 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }
      };
      unsigned int mode_5_8[8][16] = {
        { 20, 36, 34, 31, 20, 7, 41, 34, -1, 10, 41, -1, -1, -1, -1, -1 },
        { 30, 27, -1, 18, -1, 12, 20, 14, 2, 25, 15, 6, -1, -1, -1, -1 },
        { 35, -1, 41, -1, 40, -1, 39, -1, 28, -1, -1, 3, 28, -1, -1, -1 },
        { 29, -1, 0, -1, -1, 22, -1, 4, -1, 28, -1, 27, 24, 23, -1, -1 },
        { -1, 31, -1, 23, -1, 21, -1, 20, -1, 9, 12, -1, -1, 0, 13, -1 },
        { -1, 22, -1, 34, 31, -1, 14, -1, 4, -1, -1, -1, -1, -1, 22, 24 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }
      };
      unsigned int mode_1_2[8][16] = {
        { 40, -1, 38, -1, 13, -1, 5, -1, 18, -1, -1, -1, -1, -1, -1, -1 },
        { 34, -1, 35, -1, 27, -1, -1, 30, 2, 1, -1, -1, -1, -1, -1, -1 },
        { -1, 36, -1, 31, -1, 7, -1, 34, -1, 10, 41, -1, -1, -1, -1, -1 },
        { -1, 27, -1, 18, -1, 12, 20, -1, -1, -1, 15, 6, -1, -1, -1, -1 },
        { 35, -1, 41, -1, 40, -1, 39, -1, 28, -1, -1, 3, 28, -1, -1, -1 },
        { 29, -1, 0, -1, -1, 22, -1, 4, -1, 28, -1, 27, -1, 23, -1, -1 },
        { -1, 31, -1, 23, -1, 21, -1, 20, -1, -1, 12, -1, -1, 0, 13, -1 },
        { -1, 22, -1, 34, 31, -1, 14, -1, 4, -1, -1, -1, 13, -1, 22, 24 }
      };
    
      unsigned int (* mode)[8][16];
      unsigned int    n;
      unsigned int    k;
    
      if( modeNum==3 ){
        mode = &mode_13_16;
        n    = 16 * EXPANSION;
        k    = 16 * EXPANSION * 13 / 16;
      }
      else if( modeNum==2 ){
        mode = &mode_3_4;
        n    = 16 * EXPANSION;
        k    = 16 * EXPANSION * 3 / 4;
      }
      else if( modeNum==1 ){
        mode = &mode_5_8;
        n    = 16 * EXPANSION;
        k    = 16 * EXPANSION * 5 / 8;
      }
      else if( modeNum==0 ){
        mode = &mode_1_2;
        n    = 16 * EXPANSION;
        k    = 16 * EXPANSION * 1 / 2;
      }
    
      parityMatrix->initialize( n-k, n );
    
      // Set bits based on matrix definitions
      for( unsigned int rowNum=0; rowNum<(n-k)/EXPANSION; ++rowNum )
        for( unsigned int colNum=0; colNum<16; ++colNum ){
    
          for( unsigned int identityMatrix=0; identityMatrix<EXPANSION; ++identityMatrix ){
            int rotation = (*mode)[rowNum][colNum];
    
            if( rotation!=-1 ){
              unsigned int rowPos = rowNum*EXPANSION + identityMatrix;
              unsigned int colPos = colNum*EXPANSION + ((identityMatrix + rotation) % EXPANSION);
    
              parityMatrix->set( rowPos, colPos, 1 );
            }
          }
        }
    
      result->set_code( parityMatrix );
    }
    
    int main(int argc, char *argv[])
    {
      // some parameters
      int64_t maxBits     = 5000000000LL; // maximum number of bits simulated for each SNR point
      int     targetBers  = 20000;        // target number of bit errors per SNR point
      double  minBer      = 1e-10;       // BER at which to terminate simulation
      vec     EbN0db      = "0.6:0.2:5";
    
      LDPC_Code   code;
      LDPC_Parity parity;
    
      // set up space to store results
      std::vector<float> ebn0( 23, 0.0 );
    
      // Simulate each mode
      for( unsigned int modeNum=0; modeNum<4; ++modeNum ){
        std::vector< std::pair<float, float> > thisModesMeasurements;
    
        ExpandWigig( &code, &parity, modeNum );
    
        // Code settings: 40 iterations, high resolution LLR algebra
        code.set_exit_conditions( 30 );
    
        cout << code << endl;
    
        int  N      = code.get_nvar();             // number of bits per codeword
        bvec bitsin = zeros_b(N);
    
        BPSK Mod;
        vec  s = Mod.modulate_bits(bitsin);
    
        RNG_randomize();
    
        for( int j = 0; j<length(EbN0db); j++) {
          // Noise variance is N0/2 per dimension
          double       N0 = pow( 10.0, -EbN0db(j) /10.0 ) /code.get_rate();
          AWGN_Channel chan( N0 / 2 );
    
          BERC berc;  // Counters for coded and uncoded BER
    
          for( int64_t i=0; i<maxBits; i+=code.get_nvar() ){
            // Received data
            vec x = chan(s);
    
            // Demodulate
            vec softbits = Mod.demodulate_soft_bits( x, N0 );
    
            // Decode the received bits
            QLLRvec llr;
            QLLRvec llr_in = code.get_llrcalc().to_qllr(softbits);
    
            code.bp_decode( llr_in, llr );
            bvec bitsout = llr < 0;
    
            // Count the number of errors
            berc.count(bitsin, bitsout);
    
            if( berc.get_errors() > targetBers)
              break;
          }
    
          // got result for this EbN0
          thisModesMeasurements.push_back( make_pair( EbN0db(j), berc.get_errorrate() ) );
    
          cout << "Eb/N0 = " << EbN0db(j) << ": "
               << berc.get_errors()     << " errors / "
               << berc.get_total_bits() << " bits. "
               << "BER = " << berc.get_errorrate() << endl << flush;
    
          if( berc.get_errorrate() < minBer )
            break;
        }
      }
    }
    
     
    Last edit: Jason 2013-04-05
  • li li
    li li
    2013-04-05

    1, Check the definition of SNR, which may be different from the literature.

    2, Have you ever plotted the EXIT chart instead of BER? Because the error floor is caused by the intersection of CND and VND, if the shapes of them are incorrect, it may be something wrong with your LDPC decoder.

     
    Last edit: li li 2013-04-05
  • Jason
    Jason
    2013-05-02

    Thanks, Li. It looks like you're right - the errors tend to accumulate in the parity bits, so error floor excluding the parity bits is much lower (around 2*10^6, 3 orders of magnitude lower). I'll re-post the code with the modifications, just in case it's useful to someone:

    #include <iostream>
    
    #include <vector>
    #include <algorithm>
    
    #include <math.h>
    
    #include <time.h>
    
    #include <itpp/itcomm.h>
    using namespace std;
    using namespace itpp;
    
    void ExpandWigig( LDPC_Code* result, LDPC_Parity* parityMatrix, unsigned int modeNum )
    {
      const unsigned int EXPANSION = 42;
    
      unsigned int mode_13_16[8][16] = {
        { 29, 30, 0, 8, 33, 22, 17, 4, 27, 28, 20, 27, 24, 23, -1, -1 },
        { 37, 31, 18, 23, 11, 21, 6, 20, 32, 9, 12, 29, 10, 0, 13, -1 },
        { 25, 22, 4, 34, 31, 3, 14, 15, 4, 2, 14, 18, 13, 13, 22, 24 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }
      };
      unsigned int mode_3_4[8][16] = {
        { 35, 19, 41, 22, 40, 41, 39, 6, 28, 18, 17, 3, 28, -1, -1, -1 },
        { 29, 30, 0, 8, 33, 22, 17, 4, 27, 28, 20, 27, 24, 23, -1, -1 },
        { 37, 31, 18, 23, 11, 21, 6, 20, 32, 9, 12, 29, -1, 0, 13, -1 },
        { 25, 22, 4, 34, 31, 3, 14, 15, 4, -1, 14, 18, 13, 13, 22, 24 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }
      };
      unsigned int mode_5_8[8][16] = {
        { 20, 36, 34, 31, 20, 7, 41, 34, -1, 10, 41, -1, -1, -1, -1, -1 },
        { 30, 27, -1, 18, -1, 12, 20, 14, 2, 25, 15, 6, -1, -1, -1, -1 },
        { 35, -1, 41, -1, 40, -1, 39, -1, 28, -1, -1, 3, 28, -1, -1, -1 },
        { 29, -1, 0, -1, -1, 22, -1, 4, -1, 28, -1, 27, 24, 23, -1, -1 },
        { -1, 31, -1, 23, -1, 21, -1, 20, -1, 9, 12, -1, -1, 0, 13, -1 },
        { -1, 22, -1, 34, 31, -1, 14, -1, 4, -1, -1, -1, -1, -1, 22, 24 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
        { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }
      };
      unsigned int mode_1_2[8][16] = {
        { 40, -1, 38, -1, 13, -1, 5, -1, 18, -1, -1, -1, -1, -1, -1, -1 },
        { 34, -1, 35, -1, 27, -1, -1, 30, 2, 1, -1, -1, -1, -1, -1, -1 },
        { -1, 36, -1, 31, -1, 7, -1, 34, -1, 10, 41, -1, -1, -1, -1, -1 },
        { -1, 27, -1, 18, -1, 12, 20, -1, -1, -1, 15, 6, -1, -1, -1, -1 },
        { 35, -1, 41, -1, 40, -1, 39, -1, 28, -1, -1, 3, 28, -1, -1, -1 },
        { 29, -1, 0, -1, -1, 22, -1, 4, -1, 28, -1, 27, -1, 23, -1, -1 },
        { -1, 31, -1, 23, -1, 21, -1, 20, -1, -1, 12, -1, -1, 0, 13, -1 },
        { -1, 22, -1, 34, 31, -1, 14, -1, 4, -1, -1, -1, 13, -1, 22, 24 }
      };
    
      unsigned int (* mode)[8][16];
      unsigned int    n;
      unsigned int    k;
    
      if( modeNum==3 ){
        mode = &mode_13_16;
        n    = 16 * EXPANSION;
        k    = 16 * EXPANSION * 13 / 16;
      }
      else if( modeNum==2 ){
        mode = &mode_3_4;
        n    = 16 * EXPANSION;
        k    = 16 * EXPANSION * 3 / 4;
      }
      else if( modeNum==1 ){
        mode = &mode_5_8;
        n    = 16 * EXPANSION;
        k    = 16 * EXPANSION * 5 / 8;
      }
      else if( modeNum==0 ){
        mode = &mode_1_2;
        n    = 16 * EXPANSION;
        k    = 16 * EXPANSION * 1 / 2;
      }
    
      parityMatrix->initialize( n-k, n );
    
      // Set bits based on matrix definitions
      for( unsigned int rowNum=0; rowNum<(n-k)/EXPANSION; ++rowNum )
        for( unsigned int colNum=0; colNum<16; ++colNum ){
    
          for( unsigned int identityMatrix=0; identityMatrix<EXPANSION; ++identityMatrix ){
            int rotation = (*mode)[rowNum][colNum];
    
            if( rotation!=-1 ){
              unsigned int rowPos = rowNum*EXPANSION + identityMatrix;
              unsigned int colPos = colNum*EXPANSION + ((identityMatrix + rotation) % EXPANSION);
    
              parityMatrix->set( rowPos, colPos, 1 );
            }
          }
        }
    
      result->set_code( parityMatrix );
    }
    
    int main(int argc, char *argv[])
    {
      // some parameters
      int64_t maxBits     = 5000000000LL; // maximum number of bits simulated for each SNR point
      int     targetBers  = 20000;        // target number of bit errors per SNR point
      double  minBer      = 1e-10;       // BER at which to terminate simulation
      vec     EbN0db      = "0.6:0.2:5";
    
      LDPC_Code   code;
      LDPC_Parity parity;
    
      // set up space to store results
      std::vector<float> ebn0( 23, 0.0 );
    
      // Simulate each mode
      for( unsigned int modeNum=0; modeNum<4; ++modeNum ){
        std::vector< std::pair<float, float> > thisModesMeasurements;
    
        ExpandWigig( &code, &parity, modeNum );
    
        // Code settings: 40 iterations, high resolution LLR algebra
        code.set_exit_conditions( 30 );
    
        cout << code << endl;
    
        int  N      = code.get_nvar();             // number of bits per codeword
        bvec bitsin = zeros_b(N);
    
        BPSK Mod;
        vec  s = Mod.modulate_bits(bitsin);
    
        RNG_randomize();
    
        for( int j = 0; j<length(EbN0db); j++) {
          // Noise variance is N0/2 per dimension
          double       N0 = pow( 10.0, -EbN0db(j) /10.0 ) /code.get_rate();
          AWGN_Channel chan( N0 / 2 );
    
          BERC berc;  // Counters for coded and uncoded BER
          unsigned int errsNoParity = 0;
          unsigned int simmedBitsNoParity = 0;
    
          for( int64_t i=0; i<maxBits; i+=code.get_nvar() ){
            // Received data
            vec x = chan(s);
    
            // Demodulate
            vec softbits = Mod.demodulate_soft_bits( x, N0 );
    
            // Decode the received bits
            QLLRvec llr;
            QLLRvec llr_in = code.get_llrcalc().to_qllr(softbits);
    
            code.bp_decode( llr_in, llr );
            bvec bitsout = llr < 0;
    
            // Count the number of errors
            berc.count(bitsin, bitsout);
    
            for( int j=0; j<parity.get_nvar() - parity.get_ncheck(); ++j )
              if( bitsin[j]!=bitsout[j] )
                ++errsNoParity;
            simmedBitsNoParity += parity.get_nvar() - parity.get_ncheck();
    
            if( berc.get_errors() > targetBers)
              break;
          }
    
          // got result for this EbN0
          thisModesMeasurements.push_back( make_pair( EbN0db(j), berc.get_errorrate() ) );
    
          cout << "Eb/N0 = " << EbN0db(j) << ": "
               << berc.get_errors()     << " errors / "
               << berc.get_total_bits() << " bits. "
               << "BER = " << berc.get_errorrate() << endl << flush;
    
          cout << "Excluding Parity bits: " << errsNoParity << " errors / "
               << simmedBitsNoParity << " bits. BER = " << 1.0*errsNoParity/simmedBitsNoParity << endl << flush;
    
          if( berc.get_errorrate() < minBer )
            break;
        }
      }
    }