hfb.c

实现HMM算法

      hmm = ab->al_qList[q]; Nq = hmm->numStates;
      if (otprob[t][q] == NULL)/*q是模型数*/
         {
            outprob = otprob[t][q] = CreateOjsprob(&ab->abMem,Nq,S);
			/*创建在t时刻模型索引号为q的观察概率矩阵，
			outprob是array[2..Nq-1][0..S][0..M]of prob*/
            for (j=2;j<Nq;j++){
               ste=hmm->svec[j].info->pdf+1; sum = 0.0;
               outprobj = outprob[j];
			   /*取出第j个状态的概率矩阵地址，
			   outprobj是array[0..S][0..M]of prob*/
               for (s=1;s<=S;s++,ste++){
                  switch (hset->hsKind){
                  case TIEDHS:  /* SOutP deals with tied mix calculation */
                  case DISCRETEHS:
                     if (S==1) {
                        outprobj[0] = NewOtprobVec(&ab->abMem,1);
                        outprobj[0][0] = SOutP(hset,s,&ot,ste);
                     } else {
                        outprobj[s] = NewOtprobVec(&ab->abMem,1);
                        outprobj[s][0] = SOutP(hset,s,&ot,ste);
                     }
		     break;
    /* Check that PLAINHS is handled correctly this  way - efficient? */
                  case PLAINHS: 
                  case SHAREDHS: /*当前执行这里*/
			 if (S==1)
			    outprobj[0] = ShStrP(hset,ste,ot.fv[s],t,fbInfo->al_inXForm,&ab->abMem);
			    /* 计算了htkbook中的公式(1.10)ShStrP: Stream Outp calculation exploiting sharing */

		     else
		        outprobj[s] = ShStrP(hset,ste,ot.fv[s],t,fbInfo->al_inXForm,&ab->abMem);
		    break;
                  default:
                     if (S==1)
		        outprobj[0] = NULL; 
		     else
		        outprobj[s] = NULL; 
                  }
                  if (S>1)
		     sum += outprobj[s][0];
               }
               if (S>1){
                  outprobj[0][0] = sum;
                  for (s=1;s<=S;s++)
                     outprobj[s][0] = sum - outprobj[s][0];
               }
               if (trace&T_OUT && NonSkipRegion(skipstart,skipend,t)) {
                  printf(" %d. ",j); PrLog(outprobj[0][0]);
                  if (S>1){
                     printf("[ ");
                     for (s=1;s<=S;s++) PrLog(outprobj[s][0]);
                     printf("]");
                  }
               }
            }
         }
      if (trace&T_OUT && NonSkipRegion(skipstart,skipend,t)) 
         printf("\n");
   }
}

/* TraceAlphaBeta: print alpha/beta values at time t, also sum
         alpha/beta product across states at t-, t, and t+ */
static void TraceAlphaBeta(AlphaBeta *ab, int t, int startq, int endq, LogDouble pr)
{
   int i,q,Nq;
   DVector aqt,bqt;
   HLink hmm;
   double summ,sump,sum;
   
   printf("Alpha/Betas at time %d\n",t);
   summ = sump = sum = LZERO;
   for (q=startq; q<=endq; q++) {
      hmm = ab->al_qList[q]; Nq = hmm->numStates;
      printf("  Q%2d: %5s           alpha             beta\n",
             q,ab->qIds[q]->name);
      aqt = ab->alphat[q]; bqt = ab->beta[t][q];
      for (i=1;i<=Nq;i++){
         printf("                "); PrLog(aqt[i]); 
         printf("     ");        PrLog(bqt[i]);
         printf("\n");
      }
      summ = LAdd(summ,aqt[1]+bqt[1]);
      for (i=2;i<Nq;i++)
         sum = LAdd(sum,aqt[i]+bqt[i]);
      sump = LAdd(sump,aqt[Nq]+bqt[Nq]);
   }
   printf("  Sums of Products:  "); PrLog(summ-pr);
   printf("(-)   "); PrLog(sum-pr); 
   printf("    ");   PrLog(sump-pr);
   printf("(+)\n");
}
         
/* SetBeamTaper: set beam start and end points according to the minimum
           duration of the models in the current sequence */
static void SetBeamTaper(PruneInfo *p, short *qDms, int Q, int T)
{
   int q,dq,i,t;
   
   /* Set leading taper *//*当前qDms[1]和qDms[Q]=2,其余的qDms为3*/
   q=1;dq=qDms[q];i=0;/*当前dq=2*/
   for (t=1;t<=T;t++) {
      while (i==dq) {
         i=0;
         if (q<Q) q++,dq=qDms[q];
         else dq=-1;
      }
      p->qHi[t]=q;/*p->qHi[t]到当前时刻t活动的模型个数的上界，因为每个模型都有一个最小的驻留，即该模型从第1个状态到最后一个状态所经历的最小状态数，即最小帧数，p->qHi是array[1..T] of top of pruning beam,当前p->qHi[t]=13*/
      i++;
   }
   q=Q;dq=qDms[q];i=0;
   for (t=T;t>=1;t--) {
      while (i==dq) {
         i=0;
         if (q>1) q--,dq=qDms[q];
         else dq=-1;
      }
      p->qLo[t]=q;/*p->qLo是array[1..T] of bottom of pruning beam，当前 p->qLo[t]=1*/
      i++;
   }
   /*    if (trace>1) 
         for (t=1;t<=T;t++)
         printf("%d: %d to %d\n",t,p->qLo[t],p->qHi[t]);
         exit(1);*/
}


/* SetBeta: allocate and calculate beta and otprob matrices */
static LogDouble SetBeta(AlphaBeta *ab, FBInfo *fbInfo, UttInfo *utt)
{

   ParmBuf pbuf;
   int i,j,t,q,Nq,lNq=0,q_at_gMax,startq,endq;
   int S, Q, T;
   DVector bqt=NULL,bqt1,bq1t1,maxP, **beta;
   float ***outprob;
   LogDouble x,y,gMax,lMax,a,a1N=0.0;
   HLink hmm;
   PruneInfo *p;
   int skipstart, skipend;
   HMMSet *hset;
   
   hset = fbInfo->al_hset;
   skipstart = fbInfo->skipstart;
   skipend = fbInfo->skipend;
   pbuf=utt->pbuf;
   S=utt->S;
   Q=utt->Q;
   T=utt->T;
   p=ab->pInfo;
   beta=ab->beta;

   maxP = CreateDVector(&gstack, Q);   /* for calculating beam width */
  
   /* Last Column t = T，初始化beta，当t = T时的beta值 */
   p->qHi[T] = Q; endq = p->qLo[T];/*printf("\n当前p->qHi[T]=%d,endq=%d\n",p->qHi[T],endq);*/
   Setotprob(ab,fbInfo,pbuf,utt->ot,T,S,Q,endq);
   /* Setotprob: allocate and calculate otprob matrix at time t ，
   将第T帧语音概率保存到ab->otprob[T]中*/
   beta[T] = CreateBetaQ(&ab->abMem,endq,Q,Q);/*创建t = T时的beta向量*/
   /* CreateBetaQ: column of DVectors covering current beam ,
   返回v，v是DVector *v，初始时 v[q] = NULL*/
   gMax = LZERO;   q_at_gMax = 0;    /* max value of beta at time T ， LZERO为(-1.0E10)*/
   for (q=Q; q>=endq; q--){/*q是模型索引,执行的是embeded training ，
	   循环地对所有的模型在t = T时刻的beta值进行初始化，执行HTKBOOK中第129页的beta初始化公式*/
      hmm = ab->al_qList[q]; Nq = hmm->numStates;/*取出第q个模型，并把第q个模型的状态数赋给Nq*/
      bqt = beta[T][q] = NewBetaVec(&ab->abMem,Nq);/*创建t = T时第q个模型的beta向量，并赋给bqt和beta[T][q]*/
	  /* NewBetaVec: create prob vector size 1..Nq */
      bqt[Nq] = (q==Q)?0.0:beta[T][q+1][lNq]+a1N;/*当前lNq=0*//*执行HTKBOOK中第129页的beta初始化的第(1)个公式，但略有不同*/
      for (i=2;i<Nq;i++)/*此循环是完成初始化的第(2)个公式*/ 
         bqt[i] = hmm->transP[i][Nq]+bqt[Nq];/*执行HTKBOOK中第129页的beta初始化公式，因为将概率取为对数了，所以原来的乘变成了现在的加*/
      outprob = ab->otprob[T][q];/*otprob是array[1..T][1..Q][2..Nq-1][0..S][0..M]of 
	prob,T是帧数，Q是模型数（多个模型一起估计），Nq是状态数，S是流数，M是混合数*/
      x = LZERO;/*下面完成初始化的第(3)个公式*/
      for (j=2; j<Nq; j++){/*计算求和部分*/
         a = hmm->transP[1][j]; y = bqt[j];/*a为转移概率，y为观察概率b*/
         if (a>LSMALL && y > LSMALL)/*#define LSMALL (-0.5E10)   /* log values < LSMALL are set to LZERO */
            x = LAdd(x,a+outprob[j][0][0]+y);/*累加求和*/
      }
      bqt[1] = x;/*完成第（3）个公式的初始化*/
      lNq = Nq; a1N = hmm->transP[1][Nq];/*lNq保存了HMM集合中从后向前数，前一个模型的状态数，a1N保存了前一个模型中从第1个状态到最后一个状态的转移概率*/
      if (x>gMax) {
         gMax = x; q_at_gMax = q;
      }
   }/*对beta初始化完毕*/
   if (trace&T_PRU && NonSkipRegion(skipstart,skipend,T) && 
       p->pruneThresh < NOPRUNE)
      printf("%d: Beta Beam %d->%d; gMax=%f at %d\n",
             T,p->qLo[T],p->qHi[T],gMax,q_at_gMax);
   
   /* Columns T-1 -> 1 */
   for (t=T-1;t>=1;t--) {      

      gMax = LZERO;   q_at_gMax = 0;    /* max value of beta at time t */
      startq = p->qHi[t+1];
      endq = (p->qLo[t+1]==1)?1:((p->qLo[t]>=p->qLo[t+1])?p->qLo[t]:p->qLo[t+1]-1);
      while (endq>1 && ab->qDms[endq-1]==0) endq--;
      /* start end-point at top of beta beam at t+1  */
      /*  unless this is outside the beam taper.     */
      /*  + 1 to allow for state q+1[1] -> q[N]      */
      /*  + 1 for each tee model preceding endq.     */
      Setotprob(ab,fbInfo,pbuf,utt->ot,t,S,startq,endq);
      beta[t] = CreateBetaQ(&ab->abMem,endq,startq,Q);
      for (q=startq;q>=endq;q--) {
         lMax = LZERO;                 /* max value of beta in model q */
         hmm = ab->al_qList[q]; 
         Nq = hmm->numStates;
         bqt = beta[t][q] = NewBetaVec(&ab->abMem,Nq);/*t时刻，第q个模型的beta值向量*/
         bqt1 = beta[t+1][q];/*t+1时刻，第q个模型的beta值向量*/
         bq1t1 = (q==Q)?NULL:beta[t+1][q+1];/*t+1时刻，第q+1个模型的beta值向量*/
         outprob = ab->otprob[t+1][q];/*outprob为第t+1时刻时第q个模型的观察概率b*/
         bqt[Nq] = (bq1t1==NULL)?LZERO:bq1t1[1];/*bq1t1[1]为HTKBOOK第129页beta计算公式第(1)个公式第2部分的第一个加和部分,这里为该式的左端赋初值为空或bq1t1[1]*/
         if (q<startq && a1N>LSMALL)
            bqt[Nq]=LAdd(bqt[Nq],beta[t][q+1][lNq]+a1N);/*计算beta计算公式第(1)个公式第2部分*/
         for (i=Nq-1;i>1;i--){/*完成beta计算公式第(2)个公式*/
            x = hmm->transP[i][Nq] + bqt[Nq];
            if (q>=p->qLo[t+1]&&q<=p->qHi[t+1])
               for (j=2;j<Nq;j++) {
                  a = hmm->transP[i][j]; y = bqt1[j];
                  if (a>LSMALL && y>LSMALL)
                     x = LAdd(x,a+outprob[j][0][0]+y);
               }
            bqt[i] = x;
            if (x>lMax) lMax = x;
            if (x>gMax) {
               gMax = x; q_at_gMax = q;
            }
         }/*beta计算公式第(2)个公式计算完毕*/
         outprob = ab->otprob[t][q];/*下面计算beta计算公式第(3)个公式*/
         x = LZERO;
         for (j=2; j<Nq; j++){
            a = hmm->transP[1][j];
            y = bqt[j];
            if (a>LSMALL && y>LSMALL)
               x = LAdd(x,a+outprob[j][0][0]+y);
         }
         bqt[1] = x;/*beta计算公式第(3)个公式计算完毕*/
         maxP[q] = lMax;/* lMax为max value of beta in model q */
         lNq = Nq; a1N = hmm->transP[1][Nq];
      }
      while (gMax-maxP[startq] > p->pruneThresh) {/*实现剪枝*/
         beta[t][startq] = NULL;/*当全局最大beta值gMax与第startq个模型的最大beta值maxP[startq]的差距大于p->pruneThresh时将第startq个模型的beta值丢失*/
         --startq;                   /* lower startq till thresh reached */
         if (startq<1) HError(7323,"SetBeta: Beta prune failed sq < 1");
      }
      while(p->qHi[t]<startq) {        /* On taper *//*实现剪枝*/
         beta[t][startq] = NULL;
         --startq;                   /* lower startq till thresh reached */
         if (startq<1) HError(7323,"SetBeta: Beta prune failed on taper sq < 1");
      }
      p->qHi[t] = startq;
      while (gMax-maxP[endq]>p->pruneThresh){/*同上实现剪枝*/
         beta[t][endq] = NULL;
         ++endq;                   /* raise endq till thresh reached */
         if (endq>startq) {
            return(LZERO);
         }
      }
      p->qLo[t] = endq;
      if (trace&T_PRU && NonSkipRegion(skipstart,skipend,t) && 
          p->pruneThresh < NOPRUNE)
         printf("%d: Beta Beam %d->%d; gMax=%f at %d\n",
                t,p->qLo[t],p->qHi[t],gMax,q_at_gMax);
   }

   /* Finally, set total prob pr */
   utt->pr = bqt[1];/*bqt[1]保存了当q==endq且t==1时的beta值得*/

   if (utt->pr <= LSMALL) {
      return LZERO;
   }

   if (trace&T_TOP) {
      printf(" Utterance prob per frame = %e\n",utt->pr/T);
      fflush(stdout);
   }
   return utt->pr;
}

/* -------------------- Top Level of F-B Updating ---------------- */

/* CheckData: check data file consistent with HMM definition */
static void CheckData(HMMSet *hset, char *fn, BufferInfo *info, 
                      Boolean twoDataFiles) 
{
   if (info->tgtVecSize!=hset->vecSize)
      HError(7350,"CheckData: Vector size in %s[%d] is incompatible with hset [%d]",
             fn,info->tgtVecSize,hset->vecSize);
   if (!twoDataFiles){
      if (info->tgtPK != hset->pkind)
         HError(7350,"CheckData: Parameterisation in %s is incompatible with hset",
                fn);
   }
}

/* ResetStacks: Reset all stacks used by StepBack function */
static void ResetStacks(AlphaBeta *ab)
{
   ResetHeap(&ab->abMem);
}

/* StepBack: Step utterance from T to 1 calculating Beta matrix*/
static Boolean StepBack(FBInfo *fbInfo, UttInfo *utt, char * datafn)
{
   LogDouble lbeta;
   LogDouble pruneThresh;
   AlphaBeta *ab;
   PruneInfo *p;
   int qt;
   
   ResetObsCache();  
   ab = fbInfo->ab;
   pruneThresh=pruneSetting.pruneInit;/*pruneSetting.pruneInit即为我们设置的初始值，如练习中的250.0*/
   do
      {
         ResetStacks(ab);
         InitPruneStats(ab);  
         p = fbInfo->ab->pInfo;
         p->pruneThresh = pruneThresh;
         qt=CreateInsts(fbInfo,ab,utt->Q,utt->tr);
		 /* CreateInsts: create array of hmm instances for current transcription ,
		 qt中保存了总的Min no of frames to get through trans mat*/
         if (qt>utt->T) {/*当前不执行这里，utt->T是number of frames in utterance */
            if (trace&T_TOP)
               printf(" Unable to traverse %d states in %d frames\n",qt,utt->T);
            HError(-7324,"StepBack: File %s - bad data or over pruning\n",datafn);
            return FALSE;
         }
         CreateBeta(ab,utt->T);
         SetBeamTaper(p,ab->qDms,utt->Q,utt->T);/*short *ab->qDms是 array[1..Q] of minimum model duration */
         CreateOtprob(ab,utt->T);
         lbeta=SetBeta(ab,fbInfo,utt);/*lbeta中保存了通过后向算法得到的总概率P，在SetBeta(ab,fbInfo,utt)中，执行utt->pr=bqt[1],return utt->pr,bqt[1]保存了当q==endq且t==1时的beta值*/
         if (lbeta>LSMALL) break;
         pruneThresh+=pruneSetting.pruneInc;
         if (pruneThresh>pruneSetting.pruneLim || pruneSetting.pruneInc==0.0) {/*当前此判断条件为假，不执行此处*/
            if (trace&T_TOP)
               printf(" No path found in beta pass\n");
            HError(-7324,"StepBack: File %s - bad data or over pruning\n",datafn);
            return FALSE;