mp_stft.c

LastWave

    UpdateSubDict(subDictAux);
    stftAux = (STFT)subDictAux->dataContainer;
    // Update the asked sub-dictionary now
//    Printf("[H]");
    UpdateStftHarmo(stft,stftAux,dict->updateTimeIdMin,dict->updateTimeIdMax);
//     XXTerminalCursorGoBackward(3);
    break;
  case HighResStft :
    Errorf("UpdateStftSubDict : HighResStft not implemented!");
    // Update the auxiliary sub-dictionaries first
    if((subDictAux = GetStftSubDict(dict,subDict->channel,RealStft,stft->windowShape,stft->windowSize,NULL))==NULL)
      Errorf("UpdateStftSubDict (Weired) : missing auxiliary sub-dictionary RealStft for HighResStft"); 
    UpdateSubDict(subDictAux);
    stftAux = (STFT)subDictAux->dataContainer;
    if((subDictAux = GetStftSubDict(dict,subDict->channel,PhaseStft,stft->windowShape,stft->windowSize,NULL))==NULL)
      Errorf("UpdateStftSubDict (Weired) : missing auxiliary sub-dictionary PhaseStft for HighResStft"); 
    UpdateSubDict(subDictAux);
    stftAux1 = (STFT)subDictAux->dataContainer;
    if((subDictAux = GetStftSubDict(dict,subDict->channel,RealStft,stft->windowShape,stft->windowSize/scaleFactorHighRes,NULL))==NULL)
      Errorf("UpdateStftSubDict (Weired) : missing auxiliary sub-dictionary RealStft for HighResStft"); 
    UpdateSubDict(subDictAux);
    stftAux2 = (STFT)subDictAux->dataContainer;
    if((subDictAux = GetStftSubDict(dict,subDict->channel,PhaseStft,stft->windowShape,stft->windowSize/scaleFactorHighRes,NULL))==NULL)
      Errorf("UpdateStftSubDict (Weired) : missing auxiliary sub-dictionary PhaseStft for HighResStft"); 
    UpdateSubDict(subDictAux);
    stftAux3 = (STFT)subDictAux->dataContainer;
    // Update the asked sub-dictionary now
//    Printf("[h]");
    UpdateStftHighRes(stft,stftAux,stftAux1,stftAux2,stftAux3,dict->updateTimeIdMin,dict->updateTimeIdMax);
//     XXTerminalCursorGoBackward(3);
    break;
  default : Errorf("UpdateStftSubDict : type %s not implemented",StftType2Name(stft->type));
  }
  subDict->flagUpToDate = YES;
  // TODO : remove that when possible
  stft->flagUpToDate=YES;
}

char GetMaxStftSubDict(SUBDICT subDict,LISTV searchRange,LWFLOAT *pMaxValue,VALUE result) 
{
  // Search range options
  unsigned short i;
  VALUE value=NULL;
  STRVALUE str=NULL;
  LISTV  lv=NULL;
  char* rangeName=NULL;
  char *type=NULL;
  RANGE range=NULL;
  LWFLOAT min,max;
  LWFLOAT f;

  // The search limits
  char flagCausal;
  long windowSizeMin,windowSizeMax;
  long timeIdMin,timeIdMax;
  long freqIdMin,freqIdMax;

  // The location of the maximum
  LWFLOAT maxValue;
  // TODO : change to unsigned 
  int maxTimeId = 0;
  int maxFreqId = 0;
  int maxFreq0Id= 0;
  int k_freqIdMin,k_freqIdMax;
  int tmp1=0,tmp2=0;

  // Some auxiliary variables
  MOLECULE molecule=NULL;
  static MOLECULE channelMol = NULL;
  ATOM atom        =NULL;
  ATOM channelAtom =NULL;
  unsigned char channel;
  unsigned short k;
  LWFLOAT *innerProdR=NULL,*innerProdI=NULL;
  DICT dict;
  SUBDICT subDictAux = NULL;
  STFT stftComplex = NULL;
  STFT stftReal    = NULL;

  STFT stft = (STFT)(subDict->dataContainer);
  // Checking arguments
  if(subDict->flagUpToDate==NO) Errorf("GetMaxStftSubDict : subDict is out of date!");

  // Default search range
  flagCausal  = NO;
  windowSizeMin = STFT_MIN_WINDOWSIZE; windowSizeMax = STFT_MAX_WINDOWSIZE;
  timeIdMin = 0; timeIdMax = stft->signalSize-1;
  freqIdMin = 0; freqIdMax = GABOR_NYQUIST_FREQID;

  // Limiting the search range if necessary
  if(searchRange != NULL) {
    i = 0;
    while(i < searchRange->length) {
      type = GetListvNth(searchRange,i,&value,&f);
      // All search range options are &strings or &listv
      if(type==strType) {
	str = CastValue(value,STRVALUE);
	if(!strcmp(str->str,"causal")) flagCausal=YES;
	i++;continue;
      } 
      if(type!=listvType) {
	i++;continue;
      }
      lv=CastValue(value,LISTV);
      if(lv->length!=2) {
	i++;continue;
      }
      // Get the range name and the range itself
      rangeName = GetListvNthStr(lv,0);
      type = GetListvNth(lv,1,(VALUE*)&value,&f);
      if(type==rangeType) {
	range = CastValue(value,RANGE);
	if(range->size >=2) Errorf("GetMaxStftSubDict : range size should be at most 2");
	min = RangeMin(range);max=RangeMax(range);
      } else if(type == numType) {
	min = max = f;
      } else Errorf("GetMaxStftSubDict : expect a &range or &num, not a '%s'",type);
      // Treat the known fields
      if(!strcmp(rangeName,"windowSize")) {
	windowSizeMin = (long) MAX(windowSizeMin,min);windowSizeMax = (long) MIN(windowSizeMax,max);
	if(!INRANGE(windowSizeMin,stft->windowSize,windowSizeMax)) return(NO);
	i++;continue;
      }
      if(!strcmp(rangeName,"timeId")) {
	timeIdMin = (long) MAX(timeIdMin,min);timeIdMax = (long) MIN(timeIdMax,max);
	if(timeIdMin>timeIdMax) return(NO);
	i++;continue;
      }
      if(!strcmp(rangeName,"time")) {
	timeIdMin = (long) MAX(timeIdMin,Time2TimeId(stft,min));timeIdMax = (long) MIN(timeIdMax,Time2TimeId(stft,max));
	if(timeIdMin>timeIdMax) return(NO);
	i++;continue;
      }
      if(!strcmp(rangeName,"freqId")) {
	freqIdMin = (long) MAX(freqIdMin,min);freqIdMax = (long)MIN(freqIdMax,max);
	if(freqIdMin>freqIdMax) return(NO);
	i++;continue;
      }
      if(!strcmp(rangeName,"freq")) {
	freqIdMin = (long)MAX(freqIdMin,Freq2FreqId(stft,min));freqIdMax = (long)MIN(freqIdMax,Freq2FreqId(stft,max));
	if(freqIdMin>freqIdMax) return(NO);
	i++;continue;
      }
      i++;continue;
    }
  }

  // Looks for the maximum over the stft
  if(pMaxValue) *pMaxValue=0.0;
  // If the search range was empty we return NO
  if(GetStftMax(stft,flagCausal,timeIdMin,timeIdMax,freqIdMin,freqIdMax,&maxTimeId,&maxFreqId,&maxValue)==NO) return(NO);
  if(pMaxValue) *pMaxValue=maxValue;
  
  // If the result is a &listv we should fill it
  if(GetTypeValue(result)==listvType) {
    ClearListv(CastValue(result,LISTV));
    lv = TNewListv();
    AppendStr2Listv(lv,"timeId");AppendInt2Listv(lv,maxTimeId);
    AppendValue2Listv(CastValue(result,LISTV),(VALUE)lv);
    lv = TNewListv();
    AppendStr2Listv(lv,"freqId");AppendInt2Listv(lv,maxFreqId);
    AppendValue2Listv(CastValue(result,LISTV),(VALUE)lv);
    return(YES);
  } 

  // If the result is a &mol we should fill it
  dict = subDict->dict;
  if(GetTypeValue(result)==moleculeType) {
    molecule = (MOLECULE)result;
    ClearMolecule(molecule);
    SizeMolecule(molecule,stft->dimHarmo*dict->nChannels);
    // What we know so far as a 'freqId' is indeed the fundamental frequency.
    maxFreq0Id = maxFreqId;
    /* First, we fill the molecule with the first channel */
    // Next we will loop (using 'maxFreqId') on the frequency of the partials
    for(k = 0; k < stft->dimHarmo; k++) {
      atom = NewAtom();
      // First, set the time-frequency content
      CopyFieldsTFContent(stft,atom);
      atom->windowShape= stft->windowShape;
      atom->windowSize = stft->windowSize;
      atom->timeId 	= maxTimeId;
      if(stft->type==HarmoStft) {
	// Determining the box around (k+1)*freq0Id where to look for a partial
	// if it is empty we return the molecule immediately without adding a new (empty!) atom.
	if(HarmoPartialBox(stft,k,maxFreq0Id,&k_freqIdMin,&k_freqIdMax) == NO) {
	  if(k==0) Errorf("GetMaxStftSubDict : (Weird) NULL first partial");
	  DeleteAtom(atom);
	  return(YES);
	}
	// Get the best partial in the box, using the (multichannel) 'real' stft of the dictionary 
	// with the corresponding window
	if(dict->nChannels>1) channel = dict->nChannels;
	else                  channel = 0;
	if((subDictAux = GetStftSubDict(subDict->dict,channel,RealStft,stft->windowShape,stft->windowSize,NULL))==NULL)
	  Errorf("GetMaxStftSubDict (Weired) : missing auxiliary sub-dictionary RealStft for setting harmonic molecule"); 
	if(subDictAux->flagUpToDate == NO) Errorf("GetMaxStftSubDict : stftReal is out of date");
	stftReal = (STFT)subDictAux->dataContainer;
	GetStftMax(stftReal,flagCausal,maxTimeId,maxTimeId,k_freqIdMin,k_freqIdMax,&tmp1,&tmp2,&f);
	// maxFreqId is the freqId of the partial 'atom'
	maxFreqId = tmp2;
      } 
      atom->freqId = maxFreqId;

      // Then, get the complex inner product (for the first channel) using the 'complex' stft of the dictionary 
      // with the corresponding window
      channel = 0;
      if((subDictAux = GetStftSubDict(subDict->dict,channel,ComplexStft,stft->windowShape,stft->windowSize,NULL))==NULL)
	Errorf("GetMaxStftSubDict (Weired) : missing auxiliary sub-dictionary ComplexStft for setting atom"); 
      if(subDictAux->flagUpToDate == NO)
	Errorf("GetMaxStftSubDict : stftComplex is out of date");
      stftComplex = (STFT)(subDictAux->dataContainer);  
      GetStftData(stftComplex,maxTimeId,&innerProdR,&innerProdI);
      atom->coeffR= innerProdR[maxFreqId/stftComplex->fRate];
      atom->coeffI= innerProdI[maxFreqId/stftComplex->fRate];
      //	  //	  Printf("channel0 %g\n",innerProdR[maxFreqId/stftComplex->fRate]);
      // Get the phase and real coeff2
      CastAtomReal(atom);
      // Compute the HighRes coefficient if asked
      if(stft->type==HighResStft) {
	//	SetRealAtom(atom,subDict->dict,subdict->dict->signal,
	//		    flagOptimize,
	//		    NO,YES);
	Errorf("GetMaxStftSubDict : Highres atom not re-implemented");
      }
      // Now we add the atom and increase the actual dimension
      AddAtom2Molecule(molecule,atom);
    }
    /* Now, we fill the other channels if necessary */
    if(dict->nChannels > 1) {
      if(channelMol==NULL) channelMol = NewMolecule();
      SizeMolecule(channelMol,molecule->dim);
      for(channel=1; channel<dict->nChannels; channel++) {
	ClearMolecule(channelMol);
	for(k = 0; k < molecule->dim; k++) {
	  atom        = GetMoleculeAtom(molecule,0,k);
	  channelAtom = CopyAtom(atom,channelAtom);
	  /* Get the complex inner product */
	  if((subDictAux = GetStftSubDict(subDict->dict,channel,ComplexStft,stft->windowShape,stft->windowSize,NULL))==NULL)
	    Errorf("GetMaxStftSubDict (Weired) : missing auxiliary sub-dictionary ComplexStft for setting atom"); 
	  if(subDictAux->flagUpToDate == NO)
	    Errorf("GetMaxStftSubDict : stftComplex is out of date");
	  stftComplex = (STFT)(subDictAux->dataContainer);  
	  GetStftData(stftComplex,channelAtom->timeId,&innerProdR,&innerProdI);
	  channelAtom->coeffR= innerProdR[((long)(channelAtom->freqId+0.5))/stftComplex->fRate];
	  channelAtom->coeffI= innerProdI[((long)(channelAtom->freqId+0.5))/stftComplex->fRate];
	  /* Get the phase and real coeff2 */
	  CastAtomReal(channelAtom);
	  /* Compute the HighRes coefficient if asked */
	  if(stft->type==HighResStft) {
	    //	SetRealAtom(atom,subDict->dict,subdict->dict->signal,
	    //		    flagOptimize,
	    //		    NO,YES);
	    Errorf("GetMaxStftSubDict : Highres atom not re-implemented");
	  }
	  // Now we add the channelAtom to the channelMolecule */
	  AddAtom2Molecule(channelMol,channelAtom);
	}
	AddChannel2Molecule(molecule,channelMol);
      }
    }
    return(YES);
  }
  Errorf("GetMaxStftSubDict : Weired error!");
}

SubDictMethods StftMethods = {
  &GetMaxStftSubDict,
  &UpdateStftSubDict
};