kva.cpp

包括RS码的编码,硬(BM)/软(KV)译码,AWGN信道调制解调仿真. 具体采用何种编译码方案和调制解调方式可在Profile.txt文件中指定(内有详细说明). 且扩展性极好,容易向其中加入新的调

				  }
				  expb=GFMPY_V2V(expb, beta, NG);
			  }
              expa=GFMPY_V2V(expa, alpha, NG);
		  }
		  //找出j*
		  js=-1;
		  minWDeg=globalMaxDegX+(kv.K-1)*kv.MaxDY; //这是当前各g[j]最大可能的WDeg值
		  for (j=kv.MaxDY; j>=0; j--)
		  {
			  if (kv.delta[j]!=0)
			  {
				 if (kv.WD[j]<=minWDeg)
				 {
					 js=j;
					 minWDeg=kv.WD[j];
				 }
			  }
		  }
          //若J非空则调整各g[j]
		  if (js!=-1)
		  {
			 for (j=0; j<=kv.MaxDY; j++)
			 {
				 if (j!=js)
				 {
					kv.g[j].DegX=max(kv.g[j].DegX, kv.g[js].DegX);
					kv.g[j].DegY=max(kv.g[j].DegY, kv.g[js].DegY);
					maxDegX=0;  maxDegY=0;
					for (x=0; x<=kv.g[j].DegX; x++)
					for (y=0; y<=kv.g[j].DegY; y++)
					{
						kv.g[j].coef[y][x]=GFMPY_V2V(kv.g[j].coef[y][x], kv.delta[js], NG)^
							            GFMPY_V2V(kv.g[js].coef[y][x], kv.delta[j], NG);
						if (kv.g[j].coef[y][x]!=0)
						{
							if (x>=maxDegX) maxDegX=x;
							if (y>=maxDegY) maxDegY=y;
							curWDeg=x+(kv.K-1)*y;
							if (curWDeg>kv.WD[j]) kv.WD[j]=curWDeg;
						}
					}
					kv.g[j].DegX=maxDegX;
					kv.g[j].DegY=maxDegY;
				 }
			 }
			 if (++kv.g[js].DegX>globalMaxDegX) globalMaxDegX=kv.g[js].DegX;
			 _ASSERT(kv.g[js].DegX<=kv.MaxC);
			 for (x=kv.g[js].DegX; x>0; x--)
			 for (y=0; y<=kv.g[js].DegY; y++)
			 {
				 kv.g[js].coef[y][x]=GFMPY_V2V(alpha, kv.g[js].coef[y][x], NG)^
					              kv.g[js].coef[y][x-1];				 
			 }
			 for (y=0; y<=kv.g[js].DegY; y++)
			 {
				 kv.g[js].coef[y][0]=GFMPY_V2V(alpha, kv.g[js].coef[y][0], NG);
			 }
			 kv.WD[js]++;
		  }
	  }
  }
  //找出Q(X,Y)
  minWDeg=kv.MaxDX+(kv.K-1)*kv.MaxDY;//+2 ??
  js=-1;
  for (j=kv.MaxDY; j>=0; j--)
  {   
	  //curWDeg=g[j].DegX+(k-1)*g[j].DegY;
	  if (kv.WD[j]<=minWDeg)
	  {
	     minWDeg=kv.WD[j];
	     js=j;
	  }	  
  }
  return &kv.g[js];
}

/*找Q(X,Y)的所有次数小于等于D的Y-ROOT
原始算法为递归形式. 本函数中实际用栈模拟递归以提高效率
返回值为y-root的数目
*/

Polynom* YRootFactor(KVAMember& kv, const BiPolynom* Q) {
    Using_Pow2Vec_Vec2Pow_NG(kv.gf);

	int curDepth, iNZ, i, j, k, nroot;
	//curDepth:当前递归深度(也就是DFS树的搜索深度); nroot:总共找到的y-root数目
	GFE c, expa;	
    for (i=0; i<kv.K+1; i++) kv.roots[i].D=-1;
    //roots[i].D表示第i层节点是Qi(0,y)的第几个根,同时也可作为第i层节点是否扩展的标志:
	//roots[i].D=-1表示还未求Qi(0,y)的根,即当前步是由父层调用进入而非子层返回
	kv.equalroot[0]=false;

    //初始步:求出Q0(x,y)
	iNZ=-1;
    for (j=0; j<=Q->DegX; j++)
	{
	    for (i=0; i<=Q->DegY; i++)
		{
			if (Q->coef[i][j]!=0)
			{
				iNZ=j;
				break;
			}
		}
		if (iNZ!=-1) break;
	}
	if (iNZ==-1)
	{  //Q(x,y)=0, 则求y-root无意义,直接返回
	   kv.YRoots[0].D=-1;
	   return 0;
	}
	else
	{  //Q0(x,y)=<<Q(x,y)>>
	   kv.Q0[0].DegX=Q->DegX-iNZ;
	   kv.Q0[0].DegY=Q->DegY;
	   //for (i=0; i<=Q->DegY; i++) Q0[0].coef[i]=&Q->coef[i][iNZ];
	   for (i=0; i<=Q->DegY; i++)
	   {
		   for (j=0; j<=kv.Q0[0].DegX; j++) kv.Q0[0].coef[i][j]=Q->coef[i][j+iNZ];
	   }
	}

    //主递归搜索过程:
	curDepth=0;
	nroot=0;
	while (curDepth>=0)
	{
		if (kv.roots[curDepth].D==-1)
		{  //<roots[curDepth].D==-1> 表示是首次由父结点的调用进入当前过程
		   iNZ=-1;
		   for (i=0; i<=kv.Q0[curDepth].DegX; i++)
		   {
			   if (kv.Q0[curDepth].coef[0][i]!=0) {iNZ=i;  break;}
		   }
		   if (iNZ==-1 && curDepth>0 && (!kv.equalroot[curDepth]))
		   {  //<Qu(x,0)==0>, 则输出已知系数为新的y-root
			  kv.YRoots[nroot].D=curDepth-1;//当前找到的y-root的阶数
			  for (j=0; j<curDepth; j++)
			  {
				  kv.YRoots[nroot].coef[j]=kv.roots[j].coef[kv.roots[j].D];
			  }
			  nroot++;
			  kv.equalroot[curDepth]=true;
		   }

		   if (curDepth<kv.K)
		   {  //<deg[u]<D>, 求Qu(0,y)的根写入RootList, 为扩展子结点作准备
			  for (j=0; j<=Q->DegY; j++) kv.f.coef[j]=kv.Q0[curDepth].coef[j][0];
			  //!!!
			  kv.f.D=Q->DegY;
			  //若BiPolynom采用横为Y系数,纵为X系数的存放形式则可不必重排系数写入f. 待以后改进
			  kv.roots[curDepth].D=PolynomRoots(kv.roots[curDepth].coef, &kv.f);
			  _ASSERT(curDepth<=kv.K && curDepth>=0);
		   }
		   else
		   {  //边界条件2: <deg[u]!<D>, 则返回上层
			  curDepth--;
			  continue;
		   }
		}
		
		_ASSERT(curDepth<=kv.K && curDepth>=0);
		if (--kv.roots[curDepth].D<0)		   
		{  //所有根已搜索过, 返回上层
		   curDepth--;
		   continue;
		}
		kv.equalroot[curDepth+1]=(kv.roots[curDepth].coef[kv.roots[curDepth].D]!=0)? false : kv.equalroot[curDepth];
		
		//对于当前被搜索的根, 计算下一个Qu(x,y)
		kv.Q0[curDepth+1].DegX=kv.Q0[curDepth].DegX+Q->DegY-1;
		for (j=0; j<=kv.Q0[curDepth+1].DegX; j++) kv.Q0[curDepth+1].coef[0][j]=0;
		expa=1;
		for (k=0; k<=Q->DegY; k++)
		{
			for (j=0; j<kv.Q0[curDepth].DegX; j++)
				kv.Q0[curDepth+1].coef[0][j]^=GFMPY_V2V(expa, kv.Q0[curDepth].coef[k][j+1], NG);
			expa=GFMPY_V2V(expa, kv.roots[curDepth].coef[kv.roots[curDepth].D], NG);
		}
		for (i=1; i<=Q->DegY; i++)
		{
			for (j=0; j<=kv.Q0[curDepth+1].DegX; j++) kv.Q0[curDepth+1].coef[i][j]=0;
			expa=1;
     		for (k=i; k<=Q->DegY; k++)
			{
				c=GFMPY_V2V(CombMod2(k, i), expa, NG);
				for (j=kv.Q0[curDepth].DegX+i; j>=i; j--)
				{
					kv.Q0[curDepth+1].coef[i][j-1]^=GFMPY_V2V(c, kv.Q0[curDepth].coef[k][j-i], NG);
					//上面的[j-1]是考虑到Q0[curDepth+1]的首列必为全零,直接将全矩阵左移一列
					//也正因此,需将第一行的计算单独提出来
				}
				expa=GFMPY_V2V(expa, kv.roots[curDepth].coef[kv.roots[curDepth].D], NG);
			}
		}
		//Qv=<<Qu>>
		iNZ=-1;
        for (j=0; j<=kv.Q0[curDepth+1].DegX; j++)
		{
	        for (i=0; i<=Q->DegY; i++)
			{
			    if (kv.Q0[curDepth+1].coef[i][j]!=0)
				{
				   iNZ=j;
				   break;
				}
			}
		    if (iNZ!=-1) break;
		}
        if (iNZ!=-1)
		{   //最左边有iNZ个非零列,将全矩阵左移iNZ列
			//(不能用改变行首地址的方法,否则无法正确释放空间)
			for (j=0; j<=Q->DegY; j++)
			{
				for (i=0; i<=kv.Q0[curDepth+1].DegX-iNZ; i++)
					kv.Q0[curDepth+1].coef[j][i]=kv.Q0[curDepth+1].coef[j][i+iNZ];
				for (i=kv.Q0[curDepth+1].DegX-iNZ+1; i<=kv.Q0[curDepth+1].DegX; i++)
					kv.Q0[curDepth+1].coef[j][i]=0;
			}
			kv.Q0[curDepth+1].DegX-=iNZ;		    
		}

        curDepth++;
	}

    kv.YRoots[nroot].D=-1;
	return kv.YRoots;
}


bool DecodeOneWord_KVA(GFE* C, const double* reliabMat,  const RSCodeParam& rsp, void* KVARegs) {	
	KVAMember* kv=(KVAMember*)KVARegs;
	if (kv->MPAStopCdt==GivenLambda)
	{
		MultiplicityMat_Gross(*kv, reliabMat, kv->MPAThres);
	}
	else if (kv->MPAStopCdt==GivenM)
	{
		MultiplicityMat_GS(*kv, reliabMat, (int)kv->MPAThres);
	}
	else
	{
		MultiplicityMat_AA(*kv, reliabMat, (int)kv->MPAThres, kv->MPAStopCdt);		
	}
	BiPolynom* Q=Interpolation(*kv, kv->M_1D);
	Polynom* yroots=YRootFactor(*kv, Q);
	/*
#ifdef _DEBUG
  	  FILE *fp=fopen("DebugOut.txt", "w");
	  fprintf(fp, "MultiplicityMat:\n");
	  for (int i=0; i<kv->N+1; i++) {
		  for (int j=0; j<kv->N; j++) {
			  fprintf(fp, "%d ", kv->M_1D[i*kv->N+j]);
		  }
		  fprintf(fp, "\n");
	  }
	  fprintf(fp, "ReliabilityMat:\n");
	  for (i=0; i<kv->N+1; i++) {
		  for (int j=0; j<kv->N; j++) {
			  fprintf(fp, "%f ", reliabMat[i*kv->N+j]);
		  }
		  fprintf(fp, "\n");
	  }
	  fprintf(fp, "Q(x,y)(Interpolation Result):\n");
	  for (i=0; i<=Q->DegY; i++) {
		  for (int j=0; j<=Q->DegX; j++) fprintf(fp, "%d ", Q->coef[i][j]);
		  fprintf(fp, "\n");
	  }
	  fprintf(fp, "Y-Roots:\n");
	  for (i=0; yroots[i].D!=-1; i++) {
		  for (int j=0; j<=yroots[i].D; j++) fprintf(fp, "%d ", yroots[i].coef[j]);
		  fprintf(fp, "\n");
	  }
	  
	  fclose(fp);
#endif
	*/
	Using_Pow2Vec_Vec2Pow_NG(kv->gf);
	int j, i;
	GFE expa;
	double maxPr=-1.0, curPr;
	
	if (yroots[0].D==-1)
	{
		//寻找Y-ROOT失败
		for (i=0; i<kv->N; i++)
		{
			maxPr=-1.0;
			for (j=0; j<kv->Q; j++)
			{
				curPr=reliabMat[i+j*kv->N];
				if (curPr>maxPr)
				{
					maxPr=curPr;
					kv->MLCode[i]=((j==0)? 0 : Pow2Vec[j-1]);
				}
			}
		}
		for (i=0; i<kv->K; i++)
		{
			expa=1;
			C[i]=0;
			for (j=0; j<kv->N; j++)
			{
				C[i]^=GFMPY_V2V(kv->MLCode[j], expa, NG);
				if (i!=0) expa=GFMPY_V2V(expa, Pow2Vec[NG-i], NG);
			}
		}
#ifdef _DEBUG
//		faultframes++;
#endif
		return false;
	}
	
	int indMLRoot, iRoot;
	for (iRoot=0; yroots[iRoot].D!=-1; iRoot++)
	{//在所有得到的y-Roots中寻找能生成最似码字的
		curPr=1.0;//码字的可信度
		for (i=0; i<kv->N; i++)
		{//逐码元求出按当前y-Root编码得到的码字
			expa=1;
			kv->MLCode[i]=0;
			for (j=0; j<=yroots[iRoot].D; j++)
			{//按当前信息多项式求值
				kv->MLCode[i]^=GFMPY_V2V(yroots[iRoot].coef[j], expa, NG);
				expa=GFMPY_V2V(expa, Pow2Vec[i], NG);
			}

			//if (C[i]==0) curPr*=reliabMat[i];
			//else curPr*=reliabMat[(Vec2ow[C[i]]+1)*kv->N+i];
			curPr*=reliabMat[(Vec2Pow[kv->MLCode[i]]+1)*kv->N+i];
		}

		if (curPr>maxPr)
		{
		   indMLRoot=iRoot;
		   maxPr=curPr;
		}
	}
	for (j=0; j<=yroots[indMLRoot].D; j++) C[j]=yroots[indMLRoot].coef[j];
	for (j=yroots[indMLRoot].D+1; j<kv->K; j++) C[j]=0;

	return true;
}