mp_maxima.c

LastWave

  MAXIMADICT maximaDict = (MAXIMADICT)(subDict->dataContainer);
  unsigned short i;
  BOOK   book;
  unsigned long n;
  MOLECULE maxMolecule;
  MOLECULE molecule;
  //ATOM   atom;

  // Checking arguments
  if(subDict->flagUpToDate==NO) Errorf("GetMaxMaximaDictSubDict : subDict is out of date!");
  if(searchRange) Errorf("GetMaxMaximaDictSubDict : searchRange is non NULL");
  if(pMaxValue) Errorf("GetMaxMaximaDictSubDict : pMaxValue is non NULL");
  if(GetTypeValue(result)!=moleculeType) Errorf("GetMaxMaximaDictSubDict : result is not a &mol");
  
  // Looks for the maximum over the maximaDict
  maxMolecule = NULL;
  // If the search range was empty we return NO
  for(i = 0; i < maximaDict->size; i++) {
    book = maximaDict->books[i];
    for(n = 0; n < book->size; n++) {
      molecule = GetBookMolecule(book,n);
      //atom = GetMoleculeAtom(molecule,0);
      // TODO : Skip molecules which first atom is not in the search range
      //if(!INRANGE(windowSizeMin,atom->windowSize,windowSizeMax)) continue;
      //if(!INRANGE(timeIdMin,atom->timeId,timeIdMax)) continue;
      //if(!INRANGE(freqIdMin,atom->freqId,freqIdMax)) continue;
      // Case where we encounter the first atom have 
      // Case where we encounter a better molecule
      if(maxMolecule==NULL || molecule->coeff2 > maxMolecule->coeff2) {
	maxMolecule = molecule;
	continue;
      }
    }
  }
  if(maxMolecule==NULL) return(NO);

  // DEBUG : TODO : remove when this can no longer occur!
  if(isnan(maxMolecule->coeff2) || maxMolecule->coeff2 <= 0.0)  Errorf("GetMaxMaximaDictSubDict : maxMolecule coeff2=%g!",maxMolecule->coeff2);

  // Fill the result
  CopyMolecule(maxMolecule,(MOLECULE)result);
  // Memorizes the maxMolecule
  maximaDict->maxMolecule = maxMolecule;
  AddRefValue(maxMolecule);
  return(YES);
}

/******************************************************/
/*
 * 	UPDATING THE MAXIMA MOLECULES
 */
/******************************************************/
/*
 * Function to monitor the update of an atom by itself
 */
static int CompareAtoms(ATOM atomReal,ATOM atomComplex)
{
  if(atomReal->windowSize==atomComplex->windowSize &&
     (int) atomReal->timeId == (int) atomComplex->timeId &&
     (int) atomReal->freqId == (int) atomComplex->freqId)
    return(YES);
  else
    return(NO);
}


/*
 * Updating an atom, given a chosen atom. 
 */
extern void RCAtomInnerProduct(const ATOM atomR,const ATOM atomC,char flagForceNumeric,LWFLOAT *pReal,LWFLOAT *pImag);

void UpdateMaximaAtom(ATOM atom,ATOM optAtom)
{
  LWFLOAT re,im;
  
  /* Checkings */
  CheckAtomReal(atom);
  CheckAtomReal(optAtom);
  
  /*
   * Inner-product between the two atoms
   * The updated one is considered complex
   * The other is real and normalized      
   */
  RCAtomInnerProduct(optAtom,atom,NO,&re,&im);
  
  /* Case of atom ORTHOGONAL to atom : no update is needed */
  if (re == 0 && im == 0) {
    CastAtomReal(atom);
    return;
  }
  
  /* Case when we need to update */
  
  /* Update atom complex inner-product */
  atom->coeffR -= sqrt(optAtom->coeff2)*re;
  atom->coeffI -= sqrt(optAtom->coeff2)*im;
  
  /* Keeping the extreme frequencies REAL */
  if(atom->freqId == 0 || atom->freqId == GABOR_NYQUIST_FREQID)
    atom->coeffI = 0.0;
  
  /*
   * Update atom phase,coeff2,... 
   */
  CastAtomReal(atom);
  
  /* WHEN WE UPDATE AN ATOM BY ITSELF, IT SHOULD BE ALWAYS DELETED */
  if(CompareAtoms(optAtom,atom)) {
    atom->coeffR	= atom->coeffI 	= 0.0;
    atom->coeff2	= 0.0;
  }
}

/*
 * Updating a molecule given a chosen optimized molecule to remove
 */
void UpdateMaximaMolecule(MOLECULE molecule,MOLECULE optMolecule)
{
  ATOM atom,optAtom;
  unsigned short kUpdate,k;
  unsigned char channel;
  CheckMoleculeNotEmpty(molecule);
  CheckMoleculeNotEmpty(optMolecule);
  if(molecule->nChannels!=optMolecule->nChannels)
    Errorf("UpdateMaximaMolecule : bad nChannels correspondance");

  molecule->coeff2 = 0.0;
  
  /* Loop on the atoms to update */
  for(channel = 0; channel < molecule->nChannels; channel ++) {
    for(kUpdate = 0; kUpdate < molecule->dim; kUpdate++) {
      atom = GetMoleculeAtom(molecule,channel,kUpdate);
      /* Loop on the chosen atoms */
      for(k = 0; k < optMolecule->dim; k++) {
	optAtom	 = GetMoleculeAtom(optMolecule,channel,k);
	/* Updating the atoms inner-product,phase and coeff2 */
	UpdateMaximaAtom(atom,optAtom);
	/* Updating accordingly the molecule coeff2 */
	molecule->coeff2 += atom->coeff2;
      }
    }
  }
}

void UpdateMaximaDictSubDict(SUBDICT subDict)
{
  MAXIMADICT maximaDict = (MAXIMADICT)(subDict->dataContainer);
  unsigned short i;
  DICT dict = subDict->dict;
  BOOK book;
  unsigned long n;

  MOLECULE molecule;

  // If already up to date, do nothing
  if(subDict->flagUpToDate) return;

  // Some checking
  if(GetTypeValue(maximaDict)!=maximaDictType)
    Errorf("UpdateMaximaDictSubDict (Weired) : data is not maximaDict!");

  // If there is a removedMolecule (which may have been optimized) we have to update 
  // the books of local maximas well as the maxMolecule.
  if(dict->removedMolecule) {
    // DEBUG
    //    Printf("%d",maximaDict->nMaxima);
    for(i = 0; i < maximaDict->size; i++) {
      book = maximaDict->books[i];
      // The increase of 'n' is only done when the molecule is not deleted
      for(n = 0; n < book->size; ) {
	molecule = GetBookMolecule(book,n);
	// Case where this is the molecule selected through GetMax (before optimizations):
	// it has to be deleted.
	if(molecule == maximaDict->maxMolecule) {
	  DeleteMoleculeFromBook(book,n);
	  maximaDict->maxMolecule = DeleteMolecule(maximaDict->maxMolecule);
	  maximaDict->nMaxima--;
	  // Do not increase 'n'
	  continue;
	}
	UpdateMaximaMolecule(molecule,dict->removedMolecule);
	// Case when we should we delete the current molecule : don't increase 'n'
	if (molecule->coeff2 < maximaDict->threshold) {
	  DeleteMoleculeFromBook(book,n);
	  maximaDict->nMaxima--;
	  continue;
	}
	n++;
      }
    }
    // DEBUG
    //    Printf("->%d\n",maximaDict->nMaxima);
  }

  // If the maximaDict is empty, we have to (re)initialize
  if(maximaDict->nMaxima==0) {
    // First, update the necessary sub-dictionaries
    for(i = 0; i < maximaDict->size; i++) {
      // Display a clock for waiting
      switch(i%4) {
//      case 0 : Printf("\\"); Flush(); break;
//      case 1 : Printf("|"); Flush(); break;
//      case 2 : Printf("/"); Flush(); break;
//      case 3 : Printf("-"); Flush(); break;
      }
      UpdateSubDict(maximaDict->subDicts[i]);
//      Printf("\b");
    }
    // Then, compute the local maxima
    PrivateInitMaximaDict(maximaDict);
    // DEBUG
    //    Printf("Init %d\n",maximaDict->nMaxima);
    // Now all sub-dictionaries are up to date
    dict->updateTimeIdMin = dict->signalSize;
    dict->updateTimeIdMax = 0;
  }

  subDict->flagUpToDate = YES;
  // DEBUG :
  if(maximaDict->maxMolecule) Errorf("UpdateMaximaDictSubDict : maxMolecule not deleted!");
}

SubDictMethods MaximaDictMethods = {
  &GetMaxMaximaDictSubDict,
  &UpdateMaximaDictSubDict
};

/*
 * The maximaDict fields
 */
static char *maxDictBookDoc          = "{} {Gets a &listv containing all the &book of local maxima of a &maximadict.}";
static char *maxDictThresholdDoc     = "{} {Gets the threshold of a &maximadict.}";
static char *maxDictNMaximaDoc       = "{} {Gets the total number of maxima of a &maximadict.}";
void *GetMaximaDictFieldsV(MAXIMADICT maximaDict,void **arg)
{
  char *field = ARG_G_GetField(arg);
  LISTV lv;
  BOOK book;
  unsigned short i;
  /* Documentation */
  if (maximaDict == NULL) {
    if(!strcmp(field,"book"))   return(maxDictBookDoc);
    if(!strcmp(field,"thresh")) return(maxDictThresholdDoc);
    if(!strcmp(field,"nmax"))   return(maxDictNMaximaDoc);
  }
  
  if(!strcmp(field,"book")) {
    lv = TNewListv();
    for(i = 0; i < maximaDict->size;i++) {
      book=(BOOK)(maximaDict->books[i]);
      AppendValue2Listv(lv,(VALUE)book);
    }
    return(GetValueField(lv,arg));
  }
  if(!strcmp(field,"thresh")) {
    return(GetFloatField(maximaDict->threshold,arg));
  }
  if(!strcmp(field,"nmax")) {
    return(GetFloatField(maximaDict->nMaxima,arg));
  }
}

struct field fieldsMaximaDict[] = {
  "book", GetMaximaDictFieldsV, NULL, NULL, NULL,
  "thresh", GetMaximaDictFieldsV, NULL, NULL, NULL,
  "nmax", GetMaximaDictFieldsV, NULL, NULL, NULL,
  NULL, NULL, NULL, NULL, NULL
};

/*
 * The type structure for MAXIMADICT
 */

TypeStruct tsMaximaDict = {

  "{{{&maximadict} {This type is the basic type for local maxima of time-frequency dictionaries for Matching Pursuit decompositions.}}}",  /* Documentation */

  &maximaDictType,       /* The basic (unique) type name */
  NULL,     /* The GetType function */                       
  
  DeleteMaximaDict,     /* The Delete function */
  NewMaximaDict,     /* The New function */
  
  NULL,       /* The copy function */
  ClearMaximaDict,       /* The clear function */
  
  ToStrMaximaDict,       /* String conversion */
  ShortPrintMaximaDict,   /* The Print function : print the object when 'print' is called */
  PrintInfoMaximaDict,   /* The PrintInfo function : called by 'info' */

  NULL,              /* The NumExtract function : used to deal with syntax like 10a */
   
  fieldsMaximaDict,      /* The list of fields */
};

/* EOF */