hfb.c

该压缩包为最新版htk的源代码,htk是现在比较流行的语音处理软件,请有兴趣的朋友下载使用

/* ----------------------------------------------------------- */
/*                                                             */
/*                          ___                                */
/*                       |_| | |_/   SPEECH                    */
/*                       | | | | \   RECOGNITION               */
/*                       =========   SOFTWARE                  */ 
/*                                                             */
/*                                                             */
/* ----------------------------------------------------------- */
/* developed at:                                               */
/*                                                             */
/*      Speech Vision and Robotics group                       */
/*      Cambridge University Engineering Department            */
/*      http://svr-www.eng.cam.ac.uk/                          */
/*                                                             */
/*      Entropic Cambridge Research Laboratory                 */
/*      (now part of Microsoft)                                */
/*                                                             */
/* ----------------------------------------------------------- */
/*         Copyright: Microsoft Corporation                    */
/*          1995-2000 Redmond, Washington USA                  */
/*                    http://www.microsoft.com                 */
/*                                                             */
/*              2002  Cambridge University                     */
/*                    Engineering Department                   */
/*                                                             */
/*   Use of this software is governed by a License Agreement   */
/*    ** See the file License for the Conditions of Use  **    */
/*    **     This banner notice must not be removed      **    */
/*                                                             */
/* ----------------------------------------------------------- */
/*         File: HFB.c: Forward Backward routines module       */
/* ----------------------------------------------------------- */

char *hfb_version = "!HVER!HFB:   3.3 [CUED 28/04/05]";
char *hfb_vc_id = "$Id: HFB.c,v 1.3 2005/05/12 15:51:24 jal58 Exp $";

#include "HShell.h"     /* HMM ToolKit Modules */
#include "HMem.h"
#include "HMath.h"
#include "HSigP.h"
#include "HAudio.h"
#include "HWave.h"
#include "HVQ.h"
#include "HParm.h"
#include "HLabel.h"
#include "HModel.h"
#include "HTrain.h"
#include "HUtil.h"
#include "HAdapt.h"
#include "HFB.h"


/* ------------------- Trace Information ------------------------------ */
/* Trace Flags */
#define T_TOP   0001    /* Top level tracing */
#define T_OPT   0002    /* Option trace */
#define T_PRU   0004    /* pruning */
#define T_ALF   0010    /* Alpha/Beta matrices */
#define T_OCC   0020    /* Occupation Counters */
#define T_TRA   0040    /* Transition Counters */
#define T_MIX   0100    /* Mixture Weights */
#define T_OUT   0200    /* Output Probabilities */
#define T_UPD   0400    /* Model updates */
#define T_TMX  01000    /* Tied Mixture Usage */
#define T_TIM  02000    /* Time elapsed in FBFile */

static int trace         =  0;
static int skipstartInit = -1;
static int skipendInit   = -1;
static Boolean alCompLevel = FALSE;   /* align model at component level */

static ConfParam *cParm[MAXGLOBS];      /* config parameters */
static int nParm = 0;

static struct { 

   LogDouble pruneInit;      /* pruning threshold initially */
   LogDouble pruneInc;       /* pruning threshold increment */
   LogDouble pruneLim;       /* pruning threshold limit */
   float minFrwdP;           /* mix prune threshold */

} pruneSetting = { NOPRUNE, 0.0, NOPRUNE, 10.0 };

static Boolean pde = FALSE;  /* partial distance elimination */
static Boolean sharedMix = FALSE; /* true if shared mixtures */

/* ------------------------- Min HMM Duration -------------------------- */

/* Recusively calculate topological order for transition matrix */
void FindStateOrder(HLink hmm,IntVec so,int s,int *d)
{
   int p;

   so[s]=0; /* GRAY */
   for (p=1;p<hmm->numStates;p++) { 
      if (hmm->transP[p][s]>LSMALL && p!=s)
         if (so[p]<0) /* WHITE */
            FindStateOrder(hmm,so,p,d);
   }
   so[s]=++(*d); /* BLACK */
}
   

/* SetMinDurs: Set minDur field in each TrAcc */
void SetMinDurs(HMMSet *hset)
{
   HMMScanState hss;
   HLink hmm;
   TrAcc *ta;
   IntVec md,so;
   int d,nDS,i,j,k;

   NewHMMScan(hset,&hss);
   do {
      hmm = hss.hmm;
      ta = (TrAcc *)GetHook(hmm->transP);
      md = CreateIntVec(&gstack,hmm->numStates);
      so = CreateIntVec(&gstack,hmm->numStates);
      for (i=1,nDS=0;i<=hmm->numStates;i++) so[i]=md[i]=-1;
      /* Find topological order for states so that we can */
      /*  find minimum duration in single ordered pass */
      FindStateOrder(hmm,md,hmm->numStates,&nDS);
      for (i=1;i<=nDS;i++) so[md[i]]=i;
      for (i=1;i<=hmm->numStates;i++) md[i]=hmm->numStates;
      for (k=1,md[1]=0;k<=nDS;k++) {
         i=so[k];
         if (i<1 || i>hmm->numStates)  continue;
         /* Find minimum duration to state i */
         for (j=1;j<hmm->numStates;j++)
            if (hmm->transP[j][i]>LSMALL) {
               d=md[j]+((i==hmm->numStates)?0:1);
               if (d<md[i]) md[i]=d;
            }
      }
      if (nDS!=hmm->numStates) {
         char buf[8192]="";
         for (j=1;j<=hmm->numStates && strlen(buf)<4096;j++) 
            if (md[j]>=hmm->numStates)
               sprintf(buf+strlen(buf),"%d ",j);
         HError(-7332,"SetMinDurs: HMM-%s with %d/%d unreachable states ( %s)",
                HMMPhysName(hset,hmm),hmm->numStates-nDS,hmm->numStates,buf);
      }
      if (md[hmm->numStates]<0 || md[hmm->numStates]>=hmm->numStates) {
         /* Should really be an error */
         HError(-7333,"SetMinDurs: Transition matrix with discontinuity");
         ta->minDur = (hmm->transP[1][hmm->numStates]>LSMALL ? 
                       0 : 1 /*hmm->numStates-2*/ ); /* Under estimate */
      }
      else
         ta->minDur = md[hmm->numStates];
      FreeIntVec(&gstack,so); FreeIntVec(&gstack,md);
   } while (GoNextHMM(&hss));
   EndHMMScan(&hss);
}

/* ----------------------------------------------------------------------- */

/* CreateTrAcc: create an accumulator for transition counts */
static TrAcc *CreateTrAcc(MemHeap *x, int numStates)
{
   TrAcc *ta;
  
   ta = (TrAcc *) New(x,sizeof(TrAcc));
   ta->tran = CreateMatrix(x,numStates,numStates);
   ZeroMatrix(ta->tran);
   ta->occ = CreateVector(x,numStates);
   ZeroVector(ta->occ);
  
   return ta;
}

/* CreateWtAcc: create an accumulator for mixture weights */
static WtAcc *CreateWtAcc(MemHeap *x, int nMix)
{
   WtAcc *wa;
   
   wa = (WtAcc *) New(x,sizeof(WtAcc));
   wa->c = CreateVector(x,nMix);
   ZeroVector(wa->c);
   wa->occ = 0.0;
   wa->time = -1; wa->prob = NULL;
   return wa;
}

/* AttachTrAccs: attach transition accumulators to hset */
static void AttachWtTrAccs(HMMSet *hset, MemHeap *x)
{
   HMMScanState hss;
   StreamElem *ste;
   HLink hmm;
  
   NewHMMScan(hset,&hss);
   do {
      hmm = hss.hmm;
      hmm->hook = (void *)0;  /* used as numEg counter */
      if (!IsSeenV(hmm->transP)) {
         SetHook(hmm->transP, CreateTrAcc(x,hmm->numStates));
         TouchV(hmm->transP);       
      }
      while (GoNextState(&hss,TRUE)) {
         while (GoNextStream(&hss,TRUE)) {
            ste = hss.ste;
            ste->hook = CreateWtAcc(x,hss.M);
         }
      }
   } while (GoNextHMM(&hss));
   EndHMMScan(&hss);
}

/* -------------------------- Initialisation ----------------------- */
void InitFB(void)
{
   int m;
   int i;
   double d;
   Boolean b;

   Register(hfb_version,hfb_vc_id);

   for (m = 0; m < 2; ++m) {
      nParm = GetConfig(m==0 ? "HFWDBKWD" : "HFB", TRUE, cParm, MAXGLOBS);
      if (nParm>0){
         if (GetConfInt(cParm,nParm,"TRACE",&i)) trace = i;
         if (GetConfInt(cParm,nParm,"HSKIPSTART",&i)) skipstartInit = i;
         if (GetConfInt(cParm,nParm,"HSKIPEND",&i)) skipendInit = i;
         if (GetConfFlt(cParm,nParm,"PRUNEINIT", &d)) pruneSetting.pruneInit = d;
         if (GetConfFlt(cParm,nParm,"PRUNEINC", &d)) pruneSetting.pruneInc = d;
         if (GetConfFlt(cParm,nParm,"PRUNELIM", &d)) pruneSetting.pruneLim = d;
         if (GetConfFlt(cParm,nParm,"MINFORPROB", &d)) pruneSetting.minFrwdP = d;
         if (GetConfBool(cParm,nParm,"ALIGNCOMPLEVEL",&b)) alCompLevel = b;
         if (GetConfBool(cParm,nParm,"PDE",&b)) pde = b;
      }
   }
}

/* Allow tools to enable top-level tracing in HFB. Only here for historical reasons */
void SetTraceFB(void)
{
   trace |= T_TOP;
}


/* Initialise the forward backward memory stacks and make initialisations  */
void InitialiseForBack(FBInfo *fbInfo, MemHeap *x, HMMSet *hset, UPDSet uset, 
                       LogDouble pruneInit, LogDouble pruneInc, 
                       LogDouble pruneLim, float minFrwdP)
{
   int s;
   AlphaBeta *ab;
  
   fbInfo->uFlags = uset;
   fbInfo->up_hset = fbInfo->al_hset = hset;
   fbInfo->twoModels = FALSE;
   fbInfo->hsKind = hset->hsKind;
   /* Accumulators attached using AttachAccs() in HERest are overwritten
      by the following line. This is ugly and needs to be sorted out. Note:
      this function is called by HERest and HVite */
   AttachWtTrAccs(hset, x);
   SetMinDurs(hset);
   fbInfo->maxM = MaxMixInSet(hset);
   fbInfo->skipstart = skipstartInit;
   fbInfo->skipend   = skipendInit;
   for (s=1;s<=hset->swidth[0];s++)
      fbInfo->maxMixInS[s] = MaxMixInSetS(hset, s);
   fbInfo->ab = (AlphaBeta *) New(x, sizeof(AlphaBeta));
   ab = fbInfo->ab;
   CreateHeap(&ab->abMem,  "AlphaBetaFB",  MSTAK, 1, 1.0, 100000, 5000000);

   if (pruneInit < NOPRUNE) {   /* cmd line takes precedence over config file */
      pruneSetting.pruneInit = pruneInit;
      pruneSetting.pruneInc  = pruneInc;
      pruneSetting.pruneLim  = pruneLim;
   }
   if (minFrwdP < NOPRUNE)
      pruneSetting.minFrwdP  = minFrwdP;

   if (pruneSetting.pruneInit < NOPRUNE) 
      if (pruneSetting.pruneInc != 0.0)
         printf("Pruning-On[%.1f %.1f %.1f]\n", pruneSetting.pruneInit, 
                pruneSetting.pruneInc, pruneSetting.pruneLim);
      else
         printf("Pruning-On[%.1f]\n", pruneSetting.pruneInit);
   else
      printf("Pruning-Off\n");
   if (hset->numSharedMix > 0)
      sharedMix = TRUE;
   if (pde) {
      if (sharedMix)
	 HError(7399,"PDE is not compatible with shared mixtures");
      printf("Partial Distance Elimination on\n");
   }
}

/* Use a different model set for alignment */
void UseAlignHMMSet(FBInfo* fbInfo, MemHeap* x, HMMSet *al_hset)
{
   int s,S;

   /* First check 2-model mode allowed for current up_hset */
   if ((fbInfo->hsKind != PLAINHS) && (fbInfo->hsKind != SHAREDHS))
      HError(7392,"Model kind not supported for fixed alignments");
   if (al_hset->hsKind != fbInfo->hsKind)
      HRError(7392,"Different kinds in alignment and update HMM sets!");
   /* check stream compatibility */
   S = al_hset->swidth[0];
   if (S != fbInfo->up_hset->swidth[0])
      HError(7392,"Different num streams in alignment and update HMM sets!");
   for (s=1; s<=S; s++) 
      if (al_hset->swidth[s] != fbInfo->up_hset->swidth[s]) 
         HError(7392,"Stream %d widths differ in alignment and update HMM sets!",s);
   /* check update flags */
   if (fbInfo->uFlags&UPTRANS) {
      HRError(7392,"Don't update transitions on a 2-model alignment"); 
      fbInfo->uFlags = fbInfo->uFlags & ~UPTRANS;
   }
   /* check input Xforms */
   if ((al_hset->xf!=NULL) && (al_hset->xf != fbInfo->up_hset->xf))
      HError(7392,"Inconsistent input transforms for align HMMSet");
    
   /* update the global alignment set features */
   fbInfo->hsKind = al_hset->hsKind;      
   fbInfo->maxM = MaxMixInSet(al_hset);
   /* dummy accs to accomodate minDir */
   AttachWtTrAccs(al_hset, x);
   SetMinDurs(al_hset);
   /* precomps */
   if ( al_hset->hsKind == SHAREDHS)
      AttachPreComps(al_hset,al_hset->hmem);
    
   fbInfo->al_hset = al_hset;  
   fbInfo->twoModels = TRUE;
}

/* Initialise the utterance memory requirements */
void InitUttInfo( UttInfo *utt, Boolean twoFiles )
{
   CreateHeap(&utt->transStack,"transStore",MSTAK, 1, 0.5, 1000,  10000);
   CreateHeap(&utt->dataStack,"dataStore",MSTAK, 1, 0.5, 1000,  10000);
   if (twoFiles)
      CreateHeap(&utt->dataStack2,"dataStore2",MSTAK, 1, 0.5, 1000,  10000);
   utt->pbuf = NULL; utt->pbuf2 = NULL; utt->tr = NULL;
}

/* InitPruneStats: initialise pruning stats */
static void InitPruneStats(AlphaBeta *ab)
{
   PruneInfo *p;

   ab->pInfo = (PruneInfo *) New(&ab->abMem, sizeof(PruneInfo));
   p = ab->pInfo;
   p->maxBeamWidth = 0;
   p->maxAlphaBeta = LZERO;
   p->minAlphaBeta = 1.0;
}


/* -------------------------- Trace Support ----------------------- */

/* CreateTraceOcc: create the array of acc occ counts */
static void CreateTraceOcc(AlphaBeta *ab, UttInfo *utt)
{
   int q;
   Vector *occa;

   printf("\n");
   ab->occa=(Vector *)New(&ab->abMem, utt->Q*sizeof(Vector));
   occa = ab->occa;
   --occa;
   for (q=1;q<=utt->Q;q++){
      occa[q] = CreateVector(&ab->abMem, ab->al_qList[q]->numStates);
      ZeroVector(occa[q]);
   }
}

/* TraceOcc: print current accumulated occ counts for all models */
static void TraceOcc(AlphaBeta *ab, UttInfo *utt, int t)
{
   int Nq, q, i;
   Vector occaq;
   HLink hmm;
   float max;

   printf("Accumulated Occ Counts at time %d\n",t);
   for (q=1; q<=utt->Q; q++){
      occaq = ab->occa[q]; hmm = ab->al_qList[q]; Nq = hmm->numStates;
      max = 0.0;        /* ignore zero vectors */
      for (i=1;i<=Nq;i++)
         if (occaq[i]>max) max = occaq[i];
      if (max>0.0) {    /* not zero so print it */
         printf("  Q%2d: %5s", q,ab->qIds[q]->name);
         for (i=1;i<=Nq;i++) printf("%7.2f",occaq[i]);
         printf("\n");
      }
   }
}

/* SetOcct: set the global occupation count for given hmm */
static void SetOcct(HLink hmm, int q, Vector occt, Vector *occa,
                    DVector aqt, DVector bqt, DVector bq1t, LogDouble pr)
{
   int i,N;
   double x;
   Vector occaq;
   
   N=hmm->numStates;
   for (i=1;i<=N;i++) {
      x = aqt[i]+bqt[i];
      if (i==1 && bq1t != NULL && hmm->transP[1][N] > LSMALL)
         x = LAdd(x,aqt[1]+bq1t[1]+hmm->transP[1][N]);
      x -= pr;
      occt[i] = (x>MINEARG) ? exp(x) : 0.0;
   }
   if (trace&T_OCC) {
      occaq = occa[q];
      for (i=1;i<=N;i++) occaq[i] += occt[i];
   }
}


/* NonSkipRegion: returns true if t is not in the skip region */
static Boolean NonSkipRegion(int skipstart, int skipend, int t)
{
   return skipstart<1 || t<skipstart || t>skipend;
}

/* PrLog: print a log value */
void PrLog(LogDouble x)
{
   if (x<LSMALL)
      printf("       LZERO");
   else
      printf("%12.5f",x);
}



/* -------------------------------------------------------------------*/

/* GetInputObs: Get input Observations for t */
void GetInputObs( UttInfo *utt, int t, HSetKind hsKind )
{

   if (utt->twoDataFiles)
      ReadAsTable(utt->pbuf2,t-1,&(utt->ot2));
   ReadAsTable(utt->pbuf,t-1,&(utt->ot));

   if (hsKind == TIEDHS)
      if (utt->twoDataFiles)
         ReadAsTable(utt->pbuf,t-1,&(utt->ot2));

}

/* --------------------------- Forward-Backward --------------------- */



/* CheckPruning: record peak alpha.beta product and position */
static void CheckPruning(AlphaBeta *ab, int t, int skipstart, int skipend)
{
   int i,q,Nq,bestq,besti,margin;
   PruneInfo *p;
   DVector aq,bq;
   HLink hmm;
   LogDouble lx,maxL;

   bestq = besti = 0;
   maxL = LZERO;
   p = ab->pInfo;
   for (q=p->qLo[t];q<=p->qHi[t];q++){
      hmm = ab->al_qList[q]; Nq = hmm->numStates;