atom_interpolate.c

LastWave

    if (coeff2C <= coeff2M || coeff2C <= coeff2P) { 
      //Warningf("TimeInterpolation : not a strict maximum %g %g %g!",coeff2M,coeff2C,coeff2P);
      return(NO);
    } 
    ParabolaInterpolateMax(coeff2M,coeff2C,coeff2P,rate,&dTimeId); 
    *pParamVal = atom->timeId + dTimeId;
    return(YES);
    /***************************************/
    /* Interpolation on the freq parameter */
    /***************************************/
  case FreqParameter :
    /* We get the value on three points to interpolate */
    if(GetNeighborsData(atom,dict,channel,RealStft,FREQ_NEIGHBOR,&rate,&coeff2M,&coeff2C,&coeff2P,NULL,NULL,NULL)==NO) {
      //Warningf("FreqInterpolation : no neighbor available");
      return(NO);
    }
    if (coeff2C <= coeff2M || coeff2C <= coeff2P) { 
      //Warningf("FreqInterpolation : not a strict maximum %g %g %g!",coeff2M,coeff2C,coeff2P);
      return(NO);
    }
    ParabolaInterpolateMax(coeff2M,coeff2C,coeff2P,rate,&dFreqId);
    *pParamVal = atom->freqId + dFreqId;
    return(YES);
    /***************************** **********/
    /* Interpolation on the chirp parameter */
    /****************************** *********/
  case ChirpParameter :
    /* TODO : enable GetNeihborsData from already chirped atom */
    if(GetNeighborsData(atom,dict,channel,ComplexStft,FREQ_NEIGHBOR,&rate,&realM,&realC,&realP,&imagM,&imagC,&imagP)==NO) {
      //Warningf("ChirpInterpolation : no neighbor available");
      return(NO);
    }
    /* Conversion to log-amplitude and phase */
    Cartesian2LogPolar(realM,imagM,&logAmpM,&phaseM);
    Cartesian2LogPolar(realC,imagC,&logAmpC,&phaseC);
    Cartesian2LogPolar(realP,imagP,&logAmpP,&phaseP);
    ParabolaInterpolate(logAmpM,logAmpC,logAmpP,rate,&aLogAmp,&b,&c,&dFreqId);
    ParabolaInterpolate(phaseM,phaseC,phaseP,rate,&aPhi,&b,&c,&kPhi);
    /* Perform estimation based on the interpolated parabolas */
    if(EstimateChirp(atom,aLogAmp,aPhi,rate,&dChirpId,&windowSize)==NO) {
      //Warningf("ChirpInterpolation : incompatible model");
      return(NO);
    }
    freqId  = atom->freqId +dFreqId;
    chirpId = atom->chirpId+dChirpId;
    /* Check whether the estimated chirplet is aliased */
    if(!INRANGE(0,freqId-chirpId*atom->windowSize/2,GABOR_NYQUIST_FREQID) 
       || !INRANGE(0,freqId+chirpId*atom->windowSize/2,GABOR_NYQUIST_FREQID)) {
      //Warningf("ChirpInterpolation : the estimated chirplet would be aliased");
      return(NO);
    }
    *pAuxParamVal = freqId;
    *pParamVal    = chirpId;
    return(YES);
  default :
    Errorf("EstimateAtomParameter (Weired) : bad parameter type %d",paramType);
  }
}


/* Performs a series of global optimizations on a molecule, using a dictionary. */
void OptimizeMolecule(MOLECULE molecule,DICT dict,LISTV optimizations)
{
  unsigned short k;
  unsigned char channel;
  ATOM atom;

  /* Optimizations options */
  unsigned short i;
  char *type=NULL;
  VALUE value=NULL;
  LWFLOAT f;
  STRVALUE str=NULL;
  LISTV  lv=NULL;
  char* rangeName=NULL;
  RANGE range=NULL;
  
  // TODO : get rid of that ?
  char borderType      = BorderPad0;
  LWFLOAT coeffR,coeffI;
  
  SIGNAL signal = NULL;
  
  char flagTouched = NO;
  unsigned int nAtoms,nAtomsAux;
  LWFLOAT timeId,freqId,chirpId;
  LWFLOAT meanTimeId,meanFreqId,meanFreq0Id,meanFreqKId,meanChirpId;
  
  if(optimizations==NULL) return;
  
  /* Now, perform the optimizations that are asked for */
  i=0;
  while(i<optimizations->length) {
    type = GetListvNth(optimizations,i,&value,&f);
    // All optimizations options so far are &strings
    if(type!=strType) 
      Errorf("OptimizeMolecule : Weired");

    str = CastValue(value,STRVALUE);
    /* 
     * Interpolation of the time : jointly across partials and channels
     */
    if(!strcmp(str->str,"time")) {
      /* First, estimate the average time over partials and channels */
      meanTimeId = 0.0;
      nAtoms = 0;
      for(k = 0; k < molecule->dim; k++) {
	channel = 0;
	//    for(channel = 0; channel < molecule->nChannels; channel++) {
	atom = GetMoleculeAtom(molecule,channel,k);
	if(EstimateAtomParameter(atom,dict,channel,TimeParameter,&timeId,NULL)) {
	  meanTimeId+=timeId;
	  nAtoms++;
	} 
      }
      /* Then, set the new time of all these atoms */
      if(nAtoms>0) {
	meanTimeId /= nAtoms;
	meanTimeId  = QuantizeParameter(meanTimeId,AtomTimeResolution);
	for(k = 0; k < molecule->dim; k++) {
	  channel = 0;
	  //    for(channel = 0; channel < molecule->nChannels; channel++) {
	  atom = GetMoleculeAtom(molecule,channel,k);
	  atom->timeId = meanTimeId;
	  flagTouched = YES;
	}
      }
      i++;continue;
    }
    /*
     * Interpolation of the frequency : jointly over channels
     * for each partial
     */
    if(!strcmp(str->str,"freq")) {
      /* The estimate is different for each partial */
      for(k = 0; k < molecule->dim; k++) {
	/* First, estimate over the channels for the current partial */
	meanFreqId = 0.0;
	nAtoms     = 0;
	channel    = 0;
	//    for(channel = 0; channel < molecule->nChannels; channel++) {
	atom = GetMoleculeAtom(molecule,channel,k);
	if(EstimateAtomParameter(atom,dict,channel,FreqParameter,&freqId,NULL)) {
	  meanFreqId+=freqId;
	  nAtoms++;
	}  
	/* Then. set the new frequency for all channels */
	if(nAtoms>0) {
	  meanFreqId /= nAtoms;
	  meanFreqId  = QuantizeParameter(meanFreqId,AtomFreqResolution);
	  channel = 0;
	  //    for(channel = 0; channel < molecule->nChannels; channel++) {
	  atom = GetMoleculeAtom(molecule,channel,k);
	  atom->freqId = meanFreqId;
	  flagTouched = YES;
	}
      }
      i++;continue;
    }	
    /*
     * Interpolation of the chirp : jointly across partials and channels
     * the frequency is estimated jointly across channels for each partial
     */
    if(!strcmp(str->str,"chirp")) {
      /* First, estimate the average chirp over partials and channels */
      meanChirpId = 0.0;
      nAtoms = 0;
      for(k = 0; k < molecule->dim; k++) {
	/* Estimate the average frequency over channels for the current partial */
	meanFreqId = 0.0;
	nAtomsAux = 0;
	channel = 0;
	//    for(channel = 0; channel < molecule->nChannels; channel++) {
	atom = GetMoleculeAtom(molecule,channel,k);
	if(EstimateAtomParameter(atom,dict,channel,ChirpParameter,&chirpId,&freqId)) {
	  meanChirpId += chirpId;
	  meanFreqId  += freqId;
	  nAtoms++;
	  nAtomsAux++;
	} 
	
	/* Then, set the new frequency of all channels of the current partial */
	if(nAtomsAux>0) {
	  meanFreqId /= nAtomsAux;
	  meanFreqId  = QuantizeParameter(meanFreqId,AtomFreqResolution);
	  channel = 0;
	  //    for(channel = 0; channel < molecule->nChannels; channel++) {
	  atom = GetMoleculeAtom(molecule,channel,k);
	  atom->freqId = meanFreqId;
	  flagTouched = YES;
	}
      }
      
      /* Eventually, set the new chirp of all channels of all partials */
      if(nAtoms>0) {
	meanChirpId /= nAtoms;
	meanChirpId  = QuantizeParameter(meanChirpId,AtomChirpResolution);
	for(k = 0; k < molecule->dim; k++) {
	  channel = 0;
	  //    for(channel = 0; channel < molecule->nChannels; channel++) {
	  atom = GetMoleculeAtom(molecule,channel,k);
	  /* The first partial */
	  if(k == 0) {
	    meanFreq0Id   = atom->freqId;
	    atom->chirpId = meanChirpId;
	  } else {
	    /* The other partials */
	    if(meanFreq0Id==0.0) {
	      Warningf("meanFreq0Id==%g",meanFreq0Id);
	      atom->chirpId = meanChirpId*(k+1);
	    } else {
	      meanFreqKId = atom->freqId;
	      atom->chirpId = meanChirpId*meanFreqKId/meanFreq0Id;
	    }
	  }
	  flagTouched = YES;
	}
      }
      i++;continue;
    } 
    
    /* 
     * Recompute all the coefficients
     */
    if(!strcmp(str->str,"recompute")) {
      flagTouched = YES;
      i++;continue;
    }

    // DEBUG : shall we do that or just continue ?
    Errorf("OptimizeMolecule : unknown optimization '%s'",str->str);  
  }

  /* Actually recompute all coefficients */
  molecule->coeff2 = 0.0;
  for(k = 0; k < molecule->dim; k++) {
    channel = 0;
    //    for(channel = 0; channel < molecule->nChannels; channel++) {
    atom = GetMoleculeAtom(molecule,channel,k);
    signal = GetChannel(dict,channel);
    SCAtomInnerProduct(signal,atom,borderType,&coeffR,&coeffI);
    atom->coeffR = coeffR;
    atom->coeffI = coeffI;
    CastAtomReal(atom);
    molecule->coeff2 += atom->coeff2;
  } 
}

// OLD STUFF
/* We increase the octave as long as it increases the coeff2 */
static void EstimateScale(ATOM atom,DICT dict,unsigned char channel)
{
  SUBDICT subDict = NULL;
  STFT    stft    = NULL;
  SIGNAL  signal  = NULL;
  // TODO : other border types
  char borderType = BorderPad0;
  ATOM copy = NULL;
  int atomTimeIdMin,atomTimeIdMax;
  /* Checking */
  CheckAtom(atom);

  signal = GetChannel(dict,channel);
  /* First we recompute the inner product */
  SCAtomInnerProduct(signal,atom,borderType,&(atom->coeffR),&(atom->coeffI));
  CastAtomReal(atom);

  copy = CopyAtom(atom,copy);
  /*
   * We try to increase the size of the atom as long as 
   * -it is possible
   * -it increases its coeff2
   */
  while(1) {
    /* Case where the window size is too big */
    if(copy->windowSize > STFT_MAX_WINDOWSIZE) break;
    ComputeWindowSupport(copy->windowSize,copy->windowShape,copy->timeId,&atomTimeIdMin,&atomTimeIdMax);
    if(atomTimeIdMin < 0 || atomTimeIdMax >= dict->signalSize) break;
    /* Case where the window size is not available */
    if((subDict = GetStftSubDict(dict,channel,ComplexStft,copy->windowShape,copy->windowSize,NULL))==NULL) break;

    stft  = (STFT)subDict->dataContainer;
    // TODO : use flagCausal ?
    if(!INRANGE(stft->firstp,copy->timeId,stft->lastp)) break;
    /* Compute the coeff2 */
    SCAtomInnerProduct(GetChannel(dict,channel),copy,borderType,&(copy->coeffR),&(copy->coeffI));
    CastAtomReal(copy);
    /* If increasing the windowSize 'increased' the coeff2, we shall try increasing again */
    if(copy->coeff2 >= atom->coeff2) {
      atom = CopyAtom(copy,atom);
    }
    copy->windowSize *= 2;
  }
}

/* Performs a series of optimizations on an atom/a molecule, using a dictionary. */
static void OptimizeAtom(ATOM atom,DICT dict,LISTV optimizations) 
{
  /* Optimizations options */
  unsigned short i;
  char *type=NULL;
  VALUE value=NULL;
  LWFLOAT f;
  STRVALUE str=NULL;
  LISTV  lv=NULL;
  char* rangeName=NULL;
  RANGE range=NULL;

  unsigned char channel;// TODO : pass as an argument of OptimizeAtom!
  // TODO : get rid of that ?
  char borderType      = BorderPad0;
  LWFLOAT coeffR,coeffI;

  SIGNAL signal = NULL;
  static ATOM copy = NULL;

  LWFLOAT timeId,freqId,chirpId;

  if(dict->nChannels > 1) Errorf("OptimizeAtom : cannot deal with multichannel dictionaries");

  if(optimizations==NULL) return;
  /* Now, perform the optimizations that are asked for */
  // NOTE THAT THIS WILL CHANGE A LOT IN THE NEAR FUTURE!
  i=0;
  while(i<optimizations->length) {
    type = GetListvNth(optimizations,i,&value,&f);
    // All optimizations options so far are &strings
    if(type==strType) {
      str = CastValue(value,STRVALUE);
      /* Interpolation in Time and Frequency if asked */
      if(!strcmp(str->str,"freqtime")) {
	// MEANINGLESS!!!
	channel = 0;
	// TODO : when result is NO, take it into account ????
	if(EstimateAtomParameter(atom,dict,channel,FreqParameter,&freqId,NULL)==NO) freqId = atom->freqId;
	if(EstimateAtomParameter(atom,dict,channel,TimeParameter,&timeId,NULL)==NO) timeId = atom->timeId;
	// TODO : move this to another task
	timeId = QuantizeParameter(timeId,AtomTimeResolution);
	freqId = QuantizeParameter(freqId,AtomFreqResolution);

	/* Set the interpolated time and frequency */
	signal = GetChannel(dict,channel);
	copy = CopyAtom(atom,copy);
	copy->freqId = freqId; 
	copy->timeId = timeId;
	/* 
	 * This ensures compliance of copy with CheckAtom in SCAtomInnerProduct
	 * if copy->freqId==0 or copy->freqId==GABOR_NYQUIST_FREQID
	 */
	copy->coeffI = 0.0; 
	/* Recompute the inner-product */
	// TODO : move from here ??
	SCAtomInnerProduct(signal,copy,borderType,&(copy->coeffR),&(copy->coeffI));
	CastAtomReal(copy);
	if (copy->coeff2 > atom->coeff2) 
	  atom = CopyAtom(copy,atom);
	else {
	  // DEBUG
	  Warningf("Bad Interpolation");
	  PrintAtom(atom,NO);
	  PrintAtom(copy,NO);
	  SCAtomInnerProduct(signal,atom,borderType,&(atom->coeffR),&(atom->coeffI));  
	  CastAtomReal(atom);
	}
	i++;continue;
      } 
      /* Interpolation in Frequency and Chirp (and scale) if asked */
      if(!strcmp(str->str,"chirp")) {
	// MEANINGLESS!!!
	channel = 0;
	// TODO : when result is NO, take it into account
	// WARNING : this modifies atom->freqId too !
	if(EstimateAtomParameter(atom,dict,channel,ChirpParameter,&chirpId,&freqId)==NO) {
	  freqId  = atom->freqId;
	  chirpId = atom->chirpId;
	}
	QuantizeParameter(freqId,AtomFreqResolution);
	QuantizeParameter(chirpId,AtomChirpResolution);
	atom->freqId  = freqId;
	atom->chirpId = chirpId;
	// TODO : separate chirp/scale
	// DEBUG
#ifdef DEBUG_SCALE
	Printf("s : %d ",atom->windowSize);
#endif
	EstimateScale(atom,dict,channel);
#ifdef DEBUG_SCALE
	// DEBUG
	Printf("-> %d\n",atom->windowSize);
#endif
	i++;continue;
      } 

      /* Recompute the coefficient if asked */
      // TODO : do that automatically exactly when needed without asking ?
      if(!strcmp(str->str,"recompute")) {
	// MEANINGLESS!!!
	channel = 0;
	signal = GetChannel(dict,channel);
	SCAtomInnerProduct(signal,atom,borderType,&coeffR,&coeffI);
	atom->coeffR = coeffR;
	atom->coeffI = coeffI;
      } 
    }
    i++;continue;
  }

  CastAtomReal(atom);
}

void OldOptimizeMolecule(MOLECULE molecule,DICT dict,LISTV optimizations) 
{
  unsigned short k;
  unsigned char channel = 0;
  ATOM atom;

  if(optimizations==NULL) return;
  // TODO : joint optimization over all harmonics (for chirp of harmonic atoms ...)
  // TODO : joint optimization over channels

  // Optimize and recompute the molecule->coeff2
  molecule->coeff2 = 0.0;
  for(k = 0; k < molecule->dim; k++) {
    atom = GetMoleculeAtom(molecule,channel,k);
    OptimizeAtom(atom,dict,optimizations);
    molecule->coeff2 += atom->coeff2;
  }
}

/* EOF */