blaumspectra.c

根据LDPC码中码子的构造方法中的PEG算法

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream.h>
#include <fstream.h>
#include <string.h>

#include "BlaumSpectra.h"

BlaumSpectra::BlaumSpectra(const BlockCodeGenerator *blockCodeGenerator, int threshWeight) {//NULL
  N=(*blockCodeGenerator).N;
  K=(*blockCodeGenerator).K;
  M=(*blockCodeGenerator).M;
  m=(*blockCodeGenerator).m;
  T=(*blockCodeGenerator).T;

  H=new int* [M];
  for(int i=0;i<M;i++)
    H[i]=new int[N];
  for(i=0;i<M;i++){
    for(int j=0;j<N;j++)
      H[i][j]=0;
    for(j=0;j<m;j++){
      if(T[i][j]>=0)
	H[i][T[i][j]]=1;
    }
  }  
  infoBitsLocation=(*blockCodeGenerator).infoBitsLocation;

  tmpCodeWord=new int[N];
  errorPos=new int[N];

  minWeight=N-K+1;tmpMinWeight=N-K+1;
  (*this).threshWeight=threshWeight; //neglecting weight larger than or equal to 50 codewords
  numCodeWordMinWeight=new int [threshWeight];
  for(i=0;i<threshWeight;i++) numCodeWordMinWeight[i]=0;
  codeWordMinWeight=new int** [threshWeight];
  
  timeCounter=0;
  //numOcurrence=0;

  /*
  previousTestPattern=new int[N];
  previousOrder=new int[N];
  for(i=0;i<N;i++) {
    previousTestPattern[i]=0;
    previousOrder[i]=0;
  }
  */
}
BlaumSpectra::~BlaumSpectra(void) {

  delete [] tmpCodeWord;
  delete [] errorPos;
  //delete [] previousTestPattern;
  //delete [] previousOrder;

  for(int i=0;i<M;i++){
    delete [] H[i];
    H[i]=NULL;
  }
  delete [] H;
  H=NULL;

  for(i=0;i<threshWeight;i++)
    for(int j=0;j<numCodeWordMinWeight[i];j++)
      delete [] codeWordMinWeight[i][j];
  for(i=0;i<threshWeight;i++)
    if (numCodeWordMinWeight[i]!=0)
      delete [] codeWordMinWeight[i];
  delete [] codeWordMinWeight;
  delete [] numCodeWordMinWeight;
}

void BlaumSpectra::initTmpDmin() {
  tmpMinWeight=N-K+1;
}
void BlaumSpectra::setSpeedupFact(int speedupFactor) {
  speedup=speedupFactor;
}

void BlaumSpectra::erasureDecoding(int *codeWord, double *codeWordLLR){
  int i, itmp, imed;
  double tmpLLR;

  for(i=0;i<N;i++) codeWordLLR[i]=fabs(codeWordLLR[i]); 
  //cout<<"starting deBlaum::codingErasures() ..."<<endl;
  if(speedup==1) 
    itmp=decodingErasures(codeWord, 10, codeWordLLR);   
  else if(speedup==2)
    itmp=decodingErasures(codeWord, 8, codeWordLLR); 
  else 
    itmp=decodingErasures(codeWord, 6, codeWordLLR); 
  //cout<<"End one deBlaum::codingErasures()"<<endl;
}

int BlaumSpectra::decodingErasures(int *codeWord, int numNull, double *codeWordLLR) {
  int i, j, k, d, t, l, m, itmp, thresh, *I, *J, *Index, *index;
  double dtmp;
  int indica=0;

  I=new int[M];
  J=new int[N];
  Index=new int[numNull+1];
  index=new int[N];
  for(i=0;i<M;i++) I[i]=i;
  for(j=0;j<N;j++) J[j]=j;
  heapSortLeftLow(codeWordLLR, J, N, 0, N-1);

  //check whether need calculation
  /******************************************************
  k=0;
  for(i=0;i<N;i++) {
    if(previousTestPattern[i]!=codeWord[i] || J[i]!=previousOrder[i])
      k=1; break;
  }
  if(k==0) { 
    cout<<"  skip this test  ";
    delete [] I;
    I=NULL;
    delete [] J;
    J=NULL;
    delete [] Index;
    Index=NULL;
    delete [] index;
    index=NULL;
    return(1);
  }
  else {
    for(i=0;i<N;i++) {
      previousTestPattern[i]=codeWord[i];
      previousOrder[i]=J[i];
    }
  }
  *************************************************************/
  thresh=0; m=0;
  for(k=0;k<M;k++){
  POINT1:
    if(H[I[k]][J[k+thresh]]==0){
      d=-1;
      for(i=k+1;i<M;i++)
	if(H[I[i]][J[k+thresh]]==1) {d=i;break;}
      if(d==-1){//column search fail
	Index[thresh]=k+thresh;
	thresh++;
	if(thresh!=numNull+1)	goto POINT1;
	else {thresh--; goto POINT2;}
      }
      else {//SWAP d row and k row in H matrix
	j=I[k];I[k]=I[d],I[d]=j;
      }
    }
    for(i=k+1;i<M;i++){
      if(H[I[i]][J[k+thresh]]==1)
	for(j=k+thresh;j<N;j++)
	  H[I[i]][J[j]]=(H[I[i]][J[j]]+H[I[k]][J[j]])%2;
    }
  }
  Index[thresh]=k+thresh;
  //cout<<"Should not come to this point----Too Small Parameter::numEras"<<endl;
 POINT2:
  //cout<<"the number of thresh="<<thresh<<"   "<<k+thresh<<endl;
  m=0;
  for(i=0;i<Index[thresh];i++) {
    d=-1;
    for(j=0;j<thresh;j++) 
      if(i==Index[j]) {d=1; break;}
    if(d==-1) {
      index[m]=i; m++;
    }
  }
  if(m!=k) cout<<"Errors in BlaumSpectra::decodingErasures "<<m<<" "<<k+thresh<<" "<<Index[thresh]<<endl;
  //cout<<"----Point 3---"<<m<<" "<<k+numNull<<" "<<Index[numNull]<<endl; 
  t=1;
  for(i=0;i<thresh;i++) t*=2;
  double *imed=new double[t];

  //cout<<"start to tranverse all possibilities..."<<endl;
  for(i=0;i<N;i++) tmpCodeWord[i]=codeWord[i];
  for(l=0;l<t;l++){//search all permutations
    j=l;
    for(i=0;i<thresh;i++){
      tmpCodeWord[J[Index[i]]]=j%2;j/=2;
    }    
    for(i=m-1;i>=0;i--){
      itmp=0;
      for(j=index[i]+1;j<N;j++)
	itmp+=H[I[i]][J[j]]*tmpCodeWord[J[j]];
      itmp%=2;
      tmpCodeWord[J[index[i]]]=itmp;
    }
    if(l==0) {
      if(getAbsSyndrome(tmpCodeWord)!=0) {indica=1; break;} //did not find a codeword
    }
    //cout<<getAbsSyndrome(tmpCodeWord)<<" ";    
    // 
    k=cntHammDist(); //and 
    //cout<<k<<"  ";
    if(k!=0 && k<tmpMinWeight) tmpMinWeight=k;
    if(k!=0 && k<minWeight) minWeight=k;
    if(k!=0) {
      if(k<threshWeight) {
	j=compare(k);
	if(j==0) continue; //old codeword
	else if(numCodeWordMinWeight[k]==0) {
	  numCodeWordMinWeight[k]++;
	  codeWordMinWeight[k]=new int *[numCodeWordMinWeight[k]];
	  for(i=0;i<numCodeWordMinWeight[k];i++)
	    codeWordMinWeight[k][i]=new int [k];
	  for(i=0;i<k;i++) 
	    codeWordMinWeight[k][0][i]=errorPos[i];
	}
	else {
	  tmpCWMinWeight=new int *[numCodeWordMinWeight[k]+1];
	  for(i=0;i<numCodeWordMinWeight[k]+1;i++)
	    tmpCWMinWeight[i]=new int[k];
	  for(i=0;i<numCodeWordMinWeight[k];i++)
	    for(j=0;j<k;j++)
	      tmpCWMinWeight[i][j]=codeWordMinWeight[k][i][j]; //preserve the old ones
	  for(j=0;j<k;j++)
	    tmpCWMinWeight[numCodeWordMinWeight[k]][j]=errorPos[j]; //add new one
	  for(i=0;i<numCodeWordMinWeight[k];i++)
	    delete [] codeWordMinWeight[k][i];
	  delete [] codeWordMinWeight[k];
	  codeWordMinWeight[k]=tmpCWMinWeight;
	  numCodeWordMinWeight[k]++;
	}
      }
    }
  }

  delete [] I;
  I=NULL;
  delete [] J;
  J=NULL;
  delete [] Index;
  Index=NULL;
  delete [] index;
  index=NULL;
  delete [] imed;	
  imed=NULL;
  return (indica);
}


int BlaumSpectra::cntHammDist() {
  int i, j, k;
  k=0;
  for(i=0;i<N;i++) {
    if(tmpCodeWord[i]!=0) { //because the transmitted codeword is all-zero
      errorPos[k]=i;
      k++;
    }
  }
  return(k);
}

int BlaumSpectra::compare(int kk) {
  int i, j, k, l;
  l=1; // default
  for(i=0;i<numCodeWordMinWeight[kk];i++) {
    k=0;
    for(j=0;j<kk;j++) {
      if(codeWordMinWeight[kk][i][j]!=errorPos[j]) {
	k=1; break;
      }
    }
    if(k==0) {//find a repeated codeword
      l=0; break;
    }
  }
  return(l);
}

void BlaumSpectra::writeToFile(char *outNameCW, char *outNameSpect) {
  ofstream out;
  out.open(outNameCW, ios::out);
  out.setf(ios::scientific);
  for(int i=1;i<threshWeight;i++) {
    out<<"Weight="<<i<<"    num of codewords="<<numCodeWordMinWeight[i]<<endl;
    for(int j=0;j<numCodeWordMinWeight[i];j++) {
      for(int k=0;k<i;k++) {
      out<<codeWordMinWeight[i][j][k]<<" ";
      }
      out<<endl;
    }
  }
  out.close();
  //strcat(outNameCW, ".spect");
  out.open(outNameSpect, ios::out);
  out.setf(ios::scientific);
  for(i=1;i<threshWeight;i++) 
    out<<numCodeWordMinWeight[i]<<"  ";
  out<<endl;
  out.close();
}

void BlaumSpectra::print() {
  cout<<"Spectra: ";
  for(int i=0;i<threshWeight;i++)
    cout<<numCodeWordMinWeight[i]<<"("<<i<<") ";
  if(minWeight>=threshWeight) cout<<" minWeight="<<minWeight;
  cout<<endl;
}

int BlaumSpectra::getAbsSyndrome(int *codeWord){
  //for every component of the syndrome vector
  int syndrome = 0;
  int sum;
  for (int i = 0; i < M; i++){
    sum = 0;
    //compute the value of this syndrome component
    for (int j = 0; j < m; j++){
      if (T[i][j] >= 0) {
	sum += codeWord[T[i][j]];
      }
    }
    sum = sum%2;
    syndrome += sum;
  }
  return syndrome;
}

//from lowest to largest
void BlaumSpectra::prsift(double *value, int *p, int k,int i,int l)
  { int m,n, w;
    double t;
    //HEAPSORT *w;
    t=value[p[i+k]]; w=p[i+k]; n=2*(i+1)-1; m=i;
    //t=(*p[i+k]).KEY; w=p[i+k]; n=2*(i+1)-1; m=i;
    while (n<=l-1)
      { if ((n<l-1)&&(value[p[n+k]]<value[p[n+k+1]]))
	  n=n+1;
	if (t<value[p[n+k]])
          { p[m+k]=p[n+k];
            m=n; n=2*(m+1)-1;
          }
        else n=l;
      }
    p[m+k]=w;
  }

void BlaumSpectra::heapSortLeftLow(double *value, int *p,int n, int k,int m)
  { int i,l;
    if (k<0) k=0;
    if (m>n-1) m=n-1;
    l=m-k+1;
    for (i=l/2-1; i>=0; i--)
      prsift(value,p,k,i,l);
    for (i=l-1; i>=1; i--)
      { int w=p[k]; p[k]=p[i+k]; p[i+k]=w;
        prsift(value,p,k,0,i);
      }
  }