mp_maxima.c

LastWave

  Free(maxCoeffs);
  return(NULL);
}

static void ClearMaxCoeffs(MAXCOEFFS maxCoeffs)
{
  if(maxCoeffs == NULL) Errorf("ClearMaxCoeffs : NULL input");
  maxCoeffs->size = 0;
  maxCoeffs->coeff2Min = 0.0;
  maxCoeffs->coeff2Max = 0.0;
  // Note that we do NOT de-allocate the array and do not initialize it
}

static void DoubleSizeMaxCoeffs(MAXCOEFFS maxCoeffs)
{
  if(maxCoeffs == NULL)  Errorf("DoubleSizeMaxCoeffs : NULL input");
  
  // Case of a first allocation
  if(maxCoeffs->sizeAlloc == 0) {
    maxCoeffs->coeff2s   = FloatAlloc(MAXCOEFFS_SIZEMIN);
    maxCoeffs->sizeAlloc = MAXCOEFFS_SIZEMIN;
  } 
  // Case of a resize
  else {
    /* Reallocation */
    maxCoeffs->coeff2s = (LWFLOAT *)Realloc(maxCoeffs->coeff2s,2*maxCoeffs->sizeAlloc*sizeof(LWFLOAT));
    /* Update the allocated size */
    maxCoeffs->sizeAlloc *= 2;
  }
}

static void AppendCoeff2(MAXCOEFFS maxCoeffs,LWFLOAT coeff2)
{
  // Case when we need more space to store the coeff2s
  if(maxCoeffs->size == maxCoeffs->sizeAlloc) DoubleSizeMaxCoeffs(maxCoeffs);
  
  /* Appending */
  maxCoeffs->coeff2s[maxCoeffs->size] = coeff2;
  maxCoeffs->size++;
  
  /* Updating the bounds : first appended coeff2 or ... */
  if (maxCoeffs->size == 1 || coeff2 < maxCoeffs->coeff2Min) 
    maxCoeffs->coeff2Min = coeff2;
  if (coeff2 > maxCoeffs->coeff2Max) 
    maxCoeffs->coeff2Max = coeff2;
}

/**********************************************************/
/* 
 *	FIND THE RIGHT THRESHOLD IN ORDER TO GET NMAXIMA
 */
/**********************************************************/

/* Quick-Sort procedure */
static int   CmpFloat(const void* pArg1, const void* pArg2)
{
  LWFLOAT arg1,arg2;
  
  arg1 = *(LWFLOAT *)pArg1;
  arg2 = *(LWFLOAT *)pArg2;
  if(arg1 > arg2)
    return(-1);
  if(arg1 < arg2)
    return(1);
  return(0);
}

// Returns the threshold such that "exactly" 'nMaximaTarget' coeff2s are above 
// or equal to it. If the total number of coeff2s is too small, the threshold is 
// simply the minimum value of the coeff2s.
static LWFLOAT QSortFindThreshold(MAXCOEFFS maxCoeffs,unsigned long nMaximaTarget)
{
  qsort(maxCoeffs->coeff2s,maxCoeffs->size,sizeof(LWFLOAT),&CmpFloat);
  return(maxCoeffs->coeff2s[MIN(nMaximaTarget,maxCoeffs->size)-1]);
}

static unsigned long CountMaximaAboveThreshold(MAXCOEFFS maxCoeffs,LWFLOAT threshold)
{
  unsigned long i;
  unsigned long n = 0;
  for(i = 0; i < maxCoeffs->size; i++) {
    if(maxCoeffs->coeff2s[i] >= threshold)  n++;
  }
  return(n);
}

// Returns a threshold such that "approximately" 'nMaximaTarget' coeff2s are above 
// or equal to it. If the total number of coeff2s is too small, the threshold is 
// simply the minimum value of the coeff2s.
// By "approximately" it is meant that 
//   nMaximaTarget/MAXCOEFF2_TOLERANCE <= nMaxima <= nMaximaTarget*MAXCOEFF2_TOLERANCE
#define MAXCOEFF2_TOLERANCE 1.2
static LWFLOAT DichotomyFindThreshold(MAXCOEFFS maxCoeffs,unsigned long nMaximaTarget)
{
  // The target range of 
  unsigned long nMaximaTargetMin = (unsigned long) (nMaximaTarget/MAXCOEFF2_TOLERANCE);
  unsigned long nMaximaTargetMax = (unsigned long) (nMaximaTarget*MAXCOEFF2_TOLERANCE);
  // The minimum threshold and the corresponding number of maxima
  LWFLOAT	thresholdMin	= 0.0;
  unsigned long nMaximaMin	= maxCoeffs->size;
  // The maximum threshold and the corresponding number of maxima
  LWFLOAT	thresholdMax	= maxCoeffs->coeff2Max;
  unsigned long nMaximaMax	= 0;
  // The 'average' threshold and the corresponding number of maxima
  LWFLOAT	thresholdAverage;
  unsigned long nMaximaAverage;
  
  // For debug display only : count the number of iterations in the dichotomy
  unsigned long nSteps = 0;
  
  // In case there are not enough maxima available
  if(nMaximaMin < nMaximaTarget) return(thresholdMin); // TODO : replace nMaximaTarget with nMaximaTargetMax
  // There are enough maxima available, we have to find the right threshold
  // We loop until one of the extremities of the range [nMaximaMax,nMaximaMin] 
  // is in the range [nMaximaTargetMin,nMaximaTargetMax]
  while((nMaximaMax<nMaximaTargetMin) && (nMaximaMin>nMaximaTargetMax)) {
    // Count maxima above the 'average' threshold
    // TODO : look at the shape of decay of maxima to change the 'AVERAGE' 
    // (example : threshold = sqrt(thresholdMin*thresholdMax) ???)
    thresholdAverage = (thresholdMin+thresholdMax)/2;
    nMaximaAverage   = CountMaximaAboveThreshold(maxCoeffs,thresholdAverage);
    // Case 1 : the 'average' threshold is fine
    if(INRANGE(nMaximaTargetMin,nMaximaAverage,nMaximaTargetMax)) 
      return(thresholdAverage); 
    // Case 2 : the 'average' threshold is too large
    if(nMaximaAverage < nMaximaTargetMin) {
      thresholdMax = thresholdAverage;
      nMaximaMax   = nMaximaAverage;
    }
    // Case 3 : the 'average' threshold is too small
    if(nMaximaAverage > nMaximaTargetMax) {
      thresholdMin = thresholdAverage;
      nMaximaMin   = nMaximaAverage;
    }
  }
  // This code should never be reached !
  Errorf("Should return a value");
  return(0.0);
}

/* 
 * Find the local maximas from a STFT that are above a threshold
 * If STFT is not up to date, an error is generated.
 * If book!=NULL, appends to it the molecules corresponding to the found maxima if .
 * If maxCoeffs!=NULL, appends to it the value of the local maximum.
 * and thresholds them to append them in a book.
 * The molecules are neither interpolated nor chirped.
 */
static void AppendMaximaFromStftSubDict(SUBDICT subDict,char flagCausal,char flagTimeMaxima,char flagFreqMaxima,
					LWFLOAT threshold,BOOK book,MAXCOEFFS maxCoeffs)
{
  STFT   stft;
  MOLECULE molecule;
  // TODO : change to unsigned long when prototype of QuantizeRangeLarge is changed
  int timeIdMin,timeIdMax;
  unsigned long timeId;
  LWFLOAT *prevCoeffs,*coeffs,*nextCoeffs;
  LWFLOAT prevCoeff2,coeff2,nextCoeff2;
  unsigned long freqId,f;
  LWFLOAT belowCoeff2,aboveCoeff2;

  // TODO : replace with simple explicit arguments to GetMaxStftSubDict ???
  static LISTV searchRange = NULL;
  static LISTV lvTime = NULL;
  static LISTV lvFreq = NULL;
  if(searchRange==NULL) {
    searchRange = NewListv();
    lvTime = NewListv();
    lvFreq = NewListv();
    AppendStr2Listv(lvTime,"timeId");AppendInt2Listv(lvTime,0);
    AppendStr2Listv(lvFreq,"freqId");AppendInt2Listv(lvFreq,0);
    AppendValue2Listv(searchRange,(VALUE)lvTime);
    AppendValue2Listv(searchRange,(VALUE)lvFreq);
  }
  // Basic checking :
  // TODO : place in calling generic function
  if(subDict->flagUpToDate==NO)
    Errorf("AppendMaximaFromStftSubDict : subDict is out of date");
  stft = (STFT)(subDict->dataContainer);
  // TODO : remove when there is a calling generic function
  if(GetTypeValue(stft)!=stftType)
    Errorf("AppendMaximaFromStftSubDict : subDict is not of type '%s' but '%s'",stftType,GetTypeValue(stft));

  if(stft->type!=RealStft && stft->type!=HighResStft && stft->type!=HarmoStft)
    Errorf("AppendMaximaFromStftSubDict : cannot treat a stft of type '%s'",StftType2Name(stft->type));

  // Initializing the book if necessary
  if(book != NULL) {
    ClearBook(book);
    CopyFieldsTFContent(stft,book);
  }
  // The time range search 
  // TODO : simplify ?
  if(flagCausal) QuantizeRangeLarge(stft->firstp,stft->lastp,stft->tRate,&timeIdMin,&timeIdMax);
  else           QuantizeRangeLarge(0,stft->signalSize,stft->tRate,&timeIdMin,&timeIdMax);
  timeIdMin = MAX(timeIdMin,0);
  timeIdMax = MIN(timeIdMax,stft->tRate*(stft->nFrames-1));
  // Loop on time
  for(timeId = timeIdMin; timeId <= timeIdMax;timeId += stft->tRate) {
    // Get the stft 'vertical' data
    prevCoeffs    = nextCoeffs    = NULL;
    GetStftData(stft,timeId,&coeffs,NULL);
    // When not at the border, get the previous and next 'vertical' slices of data
    if(timeId != timeIdMin && timeId != timeIdMax) {
      GetStftData(stft,timeId-stft->tRate,&prevCoeffs,NULL);
      GetStftData(stft,timeId+stft->tRate,&nextCoeffs,NULL);
    }
    // Loop on frequency
    for(freqId  = stft->freqIdMin; freqId <= stft->freqIdMax; freqId += stft->fRate) {
      // Get the 'center' coeff2 at location (timeId,freqId)
      // TODO : GetStftValue(timeId,freqId) ??
      f  = (freqId-stft->freqIdMin)/stft->fRate;
      coeff2 = coeffs[f];
      // Case where we look for maxima over time and we are not at the (time) border
      if(flagTimeMaxima && prevCoeffs != NULL && nextCoeffs != NULL) {
	prevCoeff2 = prevCoeffs[f];
	nextCoeff2 = nextCoeffs[f];
	// Detection ! Should be a local maximum, and above the threshold
	if (coeff2 >= threshold && prevCoeff2 <= coeff2 && coeff2 > nextCoeff2) {
	  if(maxCoeffs) AppendCoeff2(maxCoeffs,coeff2);
	  if(book) {
	    SetListvNthFloat(lvTime,1,timeId);
	    SetListvNthFloat(lvFreq,1,freqId);
	    molecule = NewMolecule();
	    GetMaxSubDict(subDict,searchRange,NULL,(VALUE)molecule);
	    AddMolecule2Book(book,molecule);
	  }
	  // We don't want to add the molecule a second time if it is also a freq maximum!
	  continue;
	}
      }
      // Case where we look for maxima over frequency and we are not at the (freq) border
      if (flagFreqMaxima && freqId != stft->freqIdMin && freqId != stft->freqIdMax) {
	belowCoeff2 = coeffs[f-1];
	aboveCoeff2 = coeffs[f+1];
	// Detection ! Should be a local maximum, and above the threshold
	if (coeff2 >= threshold && belowCoeff2 <= coeff2 && coeff2 > aboveCoeff2) {
	  if(maxCoeffs) AppendCoeff2(maxCoeffs,coeff2);
	  if(book) {
	    SetListvNthFloat(lvTime,1,timeId);
	    SetListvNthFloat(lvFreq,1,freqId);
	    molecule = NewMolecule();
	    GetMaxSubDict(subDict,searchRange,NULL,(VALUE)molecule);
	    AddMolecule2Book(book,molecule);
	  }
	  // We don't want to add the molecule a second time if it is also a ?? maximum!
	  continue;
	}
      }
    }
  }
}

// TODO : void AppendMaximaFromSubDict(SUBDICT subdict,char flagCausal,char flagMaxType,LWFLOAT threshold,BOOK book,MAXCOEFFS maxCoeffs)
// with flagMaxType = DictMaxTime|DictMaxFreq|DictMaxScale| ....

// Inits a maximaDict : clears all local maxima and finds new ones from 
// the sub-dictionaries, which should be up to date (else an error is generated).
static void PrivateInitMaximaDict(MAXIMADICT maximaDict) {
  unsigned short i;
  BOOK book;
  SUBDICT subDict;
  static MAXCOEFFS maxCoeffs = NULL;

  // Initialize maxCoeffs
  if(maxCoeffs == NULL)  maxCoeffs = NewMaxCoeffs();
  else                   ClearMaxCoeffs(maxCoeffs);

  // First loop to get the values of the local maxima
  for(i = 0; i < maximaDict->size; i++) {
    subDict = maximaDict->subDicts[i];
    // TODO : check if we want both time and freq maxima, and if flagCausal==YES
    AppendMaximaFromStftSubDict(subDict,YES,YES,YES,0.0,NULL,maxCoeffs);
  }
  // Find the right threshold
  //  maximaDict->threshold = QSortFindThreshold(maxCoeffs,maximaDict->nMaximaTarget);
  maximaDict->threshold = DichotomyFindThreshold(maxCoeffs,maximaDict->nMaximaTarget);

  // Second loop to add the maxima greater than the threshold to the book
  maximaDict->nMaxima = 0;
  for(i = 0; i < maximaDict->size; i++) {
    subDict = maximaDict->subDicts[i];
    book    = maximaDict->books[i];
    // TODO : check if we want both time and freq maxima, and if flagCausal==YES
    AppendMaximaFromStftSubDict(subDict,YES,YES,YES,maximaDict->threshold,book,NULL);
    maximaDict->nMaxima += book->size;
    // TODO : perform optimization of molecules here if we want to
  }
}



/************************************************************************/
/*
 *    Implementation of the SUBDICT methods 
 *    for MAXIMADICT sub-dictionaries
 *  TODO : help/doc on these functions
 */
/************************************************************************/
char GetMaxMaximaDictSubDict(SUBDICT subDict,LISTV searchRange,LWFLOAT *pMaxValue,VALUE result) 
{