atom.c

LastWave

      if (coeffR >= 0) {
	*pCosPhase = 1.0;
	*pSinPhase = 0.0;
	*pPhase	 = 0.0;
      } else {
	*pCosPhase = -1.0;
	*pSinPhase = 0.0;
	*pPhase	 = 0.5;
      }
    }
    return;
  } else {
    /* Cf. PhD thesis R. Gribonval */
    /* Computing energy */
    if(pEnergy) {
      *pEnergy    = innerProd2-realGG*(coeffR*coeffR-coeffI*coeffI)+imagGG*2*coeffR*coeffI;
      *pEnergy   *= 2;
      *pEnergy   /= 1-(realGG*realGG+imagGG*imagGG);
      // DEBUG 
      if(*pEnergy<0) Errorf("ComputeRealPhaseEnergy : (Weired) negative energy!");
    }
    /* Computing phase */
    if(pPhase) {
      real = (1-realGG)*coeffR + imagGG*coeffI;
      imag = (1+realGG)*coeffI + imagGG*coeffR;
    
      norm = sqrt(real*real+imag*imag);
      *pCosPhase  = real/norm;
      *pSinPhase  = imag/norm;

      *pPhase = atan2(imag,real)/(2*M_PI);
      if (*pPhase<0)
	*pPhase = 1+*pPhase;
      if(*pPhase >= 1.0)
	*pPhase = *pPhase-1.0;
    }
  }
}

void CastAtomReal(ATOM atom)
{
  LWFLOAT phase;
  CheckAtom(atom);

  SetAtomGG(atom);
  ComputeRealPhaseEnergy(atom->coeffR,atom->coeffI,atom->realGG,atom->imagGG,
			 &phase,&(atom->cosPhase),&(atom->sinPhase),&(atom->coeff2));
  // DEBUG
  if(isnan(atom->coeff2) || atom->coeff2 < 0.0) Errorf("CastAtomReal : atom coeff2=%g!",atom->coeff2);
}

/*********************************************
 *      The main function which builds atoms 
 *********************************************
 * Build either (depending whether atomI!=NULL or not)
 * a/ a COMPLEX NORMALIZED atom (coeff2,coeffR,coeffI,cosPhase,sinPhase are not taken into account).
 * b/ a REAL atom (coeff2,cosPhase,sinPhase are used) 
 * The 'dx' and 'x0' fields of the output signals are properly set.
 *
 * If 'flagAtomSizeOnly'==YES, the output signal(s) have 'windowSize' samples corresponding to the 'support' of the atom,
 *   (the first sample is zero)
 * Else, the output signal(s) have 'signalSize' samples.
 *
 * The function built is the product of
 * -an amplitude 'sqrt(coeff2)', for REAL atoms (the amplitude is one for COMPLEX ones).
 * -a normalized window g_s(t-timeId), with shape and size determined by 'windowShape' and 'windowSize'
 * -a modulation :
 *  COMPLEX : exp(i*2*PI*(freqId*(t-timeId)/GABOR_MAX_FREQID+(chirpId/2)*(t-timeId)^2/GABOR_MAX_FREQID))
 *  or
 *  REAL      cos(2*PI*(phi+freqId*(t-timeId)/GABOR_MAX_FREQID+(chirpId/2)*(t-timeId)^2/GABOR_MAX_FREQID))
 *
 * WARNING: BORDERTYPE IS NOT  PROPERLY TREATED YET !!!!!
 ***************************************************************************/

void BuildAtom(const ATOM atom,SIGNAL atomR,SIGNAL atomI,char borderType,char flagAtomSizeOnly)
{
  /* Parameters of the atom (we make a copy to get shorter notations) */
  LWFLOAT timeId;
  LWFLOAT freqId;
  LWFLOAT chirpId;
  LWFLOAT coeff2,cosPhase,sinPhase,phase;

  /* Tabulated data */
  SIGNAL expikR 	= NULL; 
  SIGNAL expikI 	= NULL;
  SIGNAL window 	= NULL;

  static SIGNAL tempSignal = NULL;
  unsigned long reconsSize;
  // TODO : change to (unsigned) long when prototype of ComputeWindowSupport is changed ?
  int   timeIdMin,timeIdMax;

  unsigned long i;
  long index;
  double phaseArgument;
  
  LWFLOAT ishift;
  
  LWFLOAT factor;
  LWFLOAT norm2;
  
  /* Checking */
  CheckAtom(atom);
  if(atomR == NULL)               Errorf("BuildAtom : NULL atomR (Weird Error)");
  if(!BorderTypeIsOK(borderType)) Errorf("BuildAtom : unknown borderType %d",borderType);
  
  /* We copy the atom parameters to get shorter notation */
  timeId  = atom->timeId;
  freqId  = atom->freqId;
  chirpId = atom->chirpId;

  /* Complex atoms are normalized, real atoms use coeff2 and cosPhase/sinPhase */
  if(atomI==NULL) {
    coeff2    = atom->coeff2;
    cosPhase  = atom->cosPhase;
    sinPhase  = atom->sinPhase;
    phase = atan2(atom->sinPhase,atom->cosPhase)/(2*M_PI);
    if (phase<0)     phase = 1+phase;
    if(phase >= 1.0) phase = phase-1.0;
  } else {
    coeff2    = 1.0;
    cosPhase  = 1.0;
    sinPhase  = 0.0;
    phase     = 0.0;
  }

  /* Determine the size of the output signal(s) and allocate it/them */
  if (flagAtomSizeOnly == YES)     reconsSize = atom->windowSize;
  else                             reconsSize = atom->signalSize;
  SizeSignal(atomR,reconsSize,YSIG);
  ZeroSig(atomR);
  if(atomI != NULL) {
    SizeSignal(atomI,reconsSize,YSIG);
    ZeroSig(atomI);
  }
  
  /* Compute the support of the atom [timeIdMin,timeIdMax] 
   * i.e. the set of indexes (in sample coordinates of the output signal) that have to be filled.
   */
  ComputeWindowSupport(atom->windowSize,atom->windowShape,timeId,&timeIdMin,&timeIdMax);
  /* Add the first point of the atom, which value is zero */
  // TODO : remove!
  timeIdMin -=1;
  
  /* Set the dx/x0 fields of the output signals */
  if(flagAtomSizeOnly == YES) {
    /* The output signals have a time offset with respect to the input signal */
    atomR->x0 = atom->x0 + timeIdMin*atom->dx;
    atomR->dx = atom->dx;
    if(atomI != NULL){
      atomI->x0   = atomR->x0;
      atomI->dx   = atomR->dx;
    }
    /* TODO : EXPLAIN !!!! */
    ishift  = timeIdMin-timeId;
    /* The whole output signal has to be filled without any border effect */
    timeIdMin	= 0;
    timeIdMax	= reconsSize-1;
  }
  else {
    /* The output signals have the same time offset as the input signal  */
    atomR->x0 = atom->x0;
    atomR->dx = atom->dx;
    if(atomI != NULL){
      atomI->x0 = atomR->x0;
      atomI->dx = atom->dx;
    }
    /* TODO : EXPLAIN !!!! */
    ishift 	= -timeId;
    /* TODO : treat the border effects */
    switch(borderType) {
    case BorderPad0 :
      timeIdMin	= MAX(0,timeIdMin);
      timeIdMax	= MIN(reconsSize-1,timeIdMax);
      break;
    default :
      Errorf("BuildAtom : border %s not treated yet",BorderType2Name(borderType));
    }
  }
  
  /*
   * Now, we start actually building the atom
   */
  /* First, the window of the atom */
  GetTabWindow(atom->windowShape,atom->windowSize,&window);
  for(i=timeIdMin; i <= timeIdMax; i++)  atomR->Y[i] = window->Y[i-timeIdMin];
  if (atomI != NULL) {
    for(i=timeIdMin;i<=timeIdMax;i++)    atomI->Y[i] = atomR->Y[i];
  }
  /* Then, multiply by the modulation, including the chirp */
  /* We can use the tabulated exponentials only when
   * - freqId is an integer [because freqId*(i+ishift) must be integer for all i]
   * - phase is either 0 or 0.5 modulo 1
   * - chirpId is an *even* integer        [because 0.5*chirpId*(i+ishift)^2 must be integer for all i]
   */
  if (freqId != (int) freqId || ((phase != 0) && (phase != 0.5)) || chirpId/2 == (int) (chirpId/2)) {
    /* Modulating without tabulated data */
    for(i=timeIdMin;i<=timeIdMax;i++) {
      /* TODO : EXPLAIN the use of ishift */
      phaseArgument = 2*M_PI*(phase+ (freqId*(i+ishift))/(GABOR_MAX_FREQID)
			      +0.5*(chirpId*(i+ishift)*(i+ishift))/(GABOR_MAX_FREQID));
      atomR->Y[i] *= cos(phaseArgument);
      if(atomI != NULL)
	atomI->Y[i] *= sin(phaseArgument);
    }
  } else {
    /* Modulating with tabulated data */
    GetTabExponentials(&expikR,&expikI);

    /* Computations are modulo 'GABOR_MAX_FREQID', with positive integers */
    ishift   = GABOR_MAX_FREQID + ((int) ishift)%(GABOR_MAX_FREQID);
    for(i=timeIdMin;i<=timeIdMax;i++)  {	
      /* TODO : EXPLAIN the use of ishift */
      // Is there a problem here with phase?
      //index = ((int) (phase*GABOR_MAX_FREQID
      //+ freqId*(i+ishift)
      //+ 0.5*chirpId*(i+ishift)*(i+ishift)))	%(GABOR_MAX_FREQID);
      index  = ((int)freqId)*(i+(int)ishift);
      index += ((int)(chirpId/2))*(i+(int)ishift)*(i+(int)ishift);
      if(phase==0.5)
	index += GABOR_MAX_FREQID/2;
      index = index%GABOR_MAX_FREQID;
      index = (index+GABOR_MAX_FREQID)%GABOR_MAX_FREQID;
      if(!INRANGE(0,index,expikR->size-1)) Errorf("BuildAtom : (Weird) internal accessing expikR[%d]",index);
      if(i>=atomR->size)      Errorf("BuildAtom : (Weird) internal atomR");
      atomR->Y[i] *= expikR->Y[index];
      if(atomI != NULL)
	atomI->Y[i] *= expikI->Y[index];
    }
  }
  
  /*
   * If the atom we built is complex, we are done
   */
  if(atomI!=NULL) return;
  else {
    /* Compute the norm of the real atom */
    norm2 = 0.0;
    for(i=timeIdMin;i<=timeIdMax;i++) norm2 += atomR->Y[i]*atomR->Y[i];  
    /* Normalize the atom and set its amplitude */
    factor = sqrt(coeff2/norm2);
    for(i=timeIdMin;i<=timeIdMax;i++) atomR->Y[i] *= factor;
  }
}


/*
 * The functions to Get the atom fields
 */

/*
 * The window fields : windowShape,windowSize,
 */
static char *windowShapeDoc = "{[= <windowShape>]} {Sets/Gets the windowShape of an atom. "
WindowShapeHelpString
"}";
static char *windowSizeDoc = "{[= <windowSize>]} {Sets/Gets the windowSize of an atom, i.e. the number of samples of its window. So far only powers of 2 are allowed.}";

void *GetWindowShapeAtomV(ATOM atom, void **arg)
{
  /* Documentation */
  if (atom == NULL) return(windowShapeDoc);
  return(GetStrField(WindowShape2Name(atom->windowShape),arg));
}

void *GetWindowSizeAtomV(ATOM atom, void **arg)
{
  /* Documentation */
  if (atom == NULL) return(windowSizeDoc);
  return(GetIntField(atom->windowSize,arg));
}

void *SetWindowShapeAtomV(ATOM atom, void **arg)
{
  static char *windowShapeName = NULL;
  char windowShape;
  /* Documentation */
  if (atom == NULL) return(windowShapeDoc);

  // Allocation (once only) of the resulting string
  if(windowShapeName == NULL) {
    windowShapeName = CharAlloc(1);
    windowShapeName[0] = '\0';
  }
  if(SetStrField(&windowShapeName,arg)==NULL) return(NULL);
  windowShape = Name2WindowShape(windowShapeName);
  // If the windowShape changed we have to recompute the 'GG' of the atom
  // to keep coherent
  if(atom->windowShape == windowShape) return(strType);
  atom->windowShape = windowShape;
  SetAtomGG(atom);
  // DEBUG TODO : remove!
  if(atom->realGG*atom->realGG+atom->imagGG*atom->imagGG>1) {
    PrintInfoValue(atom);
    Errorf("SetWindowShapeAtomV : (Weired) bad GG (%g %g)",atom->realGG,atom->imagGG);
  }
  return(strType);
}

void *SetWindowSizeAtomV(ATOM atom, void **arg)
{
  int windowSize;
  /* Documentation */
  if (atom == NULL) return(windowShapeDoc);
  // Init (for += syntax for example)
  windowSize = atom->windowSize;
  if(SetIntField(&windowSize,arg,FieldSPositive)==NULL) return(NULL);

  /* Some checkings */
  if(!INRANGE(STFT_MIN_WINDOWSIZE,windowSize,STFT_MAX_WINDOWSIZE) || windowSize > atom->signalSize) {
    SetErrorf("Cannot set atom.windowSize=%d. Should be between %d and %d and at most atom.signalSize (%d)",
	      windowSize,STFT_MIN_WINDOWSIZE,STFT_MAX_WINDOWSIZE,atom->signalSize);
    return(NULL);
  }
  // If the windowSize changed we have to recompute the 'GG' of the atom
  // to keep coherent
  if(atom->windowSize == windowSize) return(numType);
  atom->windowSize = windowSize;
  SetAtomGG(atom);
  // DEBUG TODO : remove!
  if(atom->realGG*atom->realGG+atom->imagGG*atom->imagGG>1) {
    PrintInfoValue(atom);
    Errorf("SetWindowSizeAtomV : (Weired) bad GG (%g %g)",atom->realGG,atom->imagGG);
  }
  return(numType);
}


/*
 * The dt/df fields
 */
static char *dtDoc = "{} {Gets the time spread of an atom in seconds.}";
static char *dfDoc = "{} {Gets the frequency spread of an atom in Hertz.}";
static char *supportDoc = "{} {Gets the time support {timeMin timeMax} of an atom in seconds}"; 
static char *supportIdDoc = "{} {Gets the time support {timeIdMin timeIdMax} of an atom in sample coordinates}"; 

void *GetDtDfAtomV(ATOM atom,void **arg)
{
  char *field = ARG_G_GetField(arg);
  LISTV lv;
  int timeIdMin,timeIdMax;
  /* Documentation */
  if (atom == NULL) {
    if(!strcmp(field,"dt")) return(dtDoc);
    if(!strcmp(field,"df")) return(dfDoc);
    if(!strcmp(field,"support")) return(supportDoc);
    if(!strcmp(field,"supportId")) return(supportIdDoc);
  }
  if(!strcmp(field,"dt"))  
    return(GetFloatField(0.5*(TimeId2Time(atom,atom->windowSize)-TimeId2Time(atom,0)),arg));
  if(!strcmp(field,"df"))
    return(GetFloatField(0.5*FreqId2Freq(atom,((LWFLOAT)GABOR_MAX_FREQID)/(2*M_PI*theGaussianSigma2*(LWFLOAT)atom->windowSize)),arg));
  if(!strcmp(field,"support") || !strcmp(field,"supportId")) {
    lv = TNewListv();
    ComputeWindowSupport(atom->windowSize,atom->windowShape,atom->timeId,
			 &timeIdMin,&timeIdMax);
    if(!strcmp(field,"supportId")) {
      AppendInt2Listv(lv,timeIdMin);
      AppendInt2Listv(lv,timeIdMax);
    } else {
      AppendFloat2Listv(lv,TimeId2Time(atom,timeIdMin));
      AppendFloat2Listv(lv,TimeId2Time(atom,timeIdMax));
    }
    return(GetValueField(lv,arg));
  }
}

/*
 * The time(Id) fields
 */
static char *timeIdDoc = "{[= <timeId>]} {Sets/Gets the time center of an atom in sample coordinates, i.e. an index 0 <= timeId < atom.signalSize.}";
static char *timeDoc = "{[= <time>]} {Sets/Gets the time center of an atom in real coordinates, i.e. the real time in seconds.}";

void *GetTimeAtomV(ATOM atom,void **arg)
{
  char *field = ARG_G_GetField(arg);
 
  /* Documentation */
  if (atom == NULL) {
    if(!strcmp(field,"time")) return(timeDoc);
    if(!strcmp(field,"timeId")) return(timeIdDoc);
  }
  if(!strcmp(field,"time"))  
    return(GetFloatField(TimeId2Time(atom,atom->timeId),arg));
  if(!strcmp(field,"timeId"))
    return(GetFloatField(atom->timeId,arg));
}

void *SetTimeAtomV(ATOM atom,void **arg)
{
  char *field = ARG_G_GetField(arg);
  char flagId = NO;

  LWFLOAT time,timeId;
  /* Documentation */
  if (atom == NULL) {
    if(!strcmp(field,"time")) return(timeDoc);
    if(!strcmp(field,"timeId")) return(timeIdDoc);
  }
  if(!strcmp(field,"timeId"))    flagId = YES;
  else if(!strcmp(field,"time")) flagId = NO;
  else Errorf("SetTimeAtomV : Weird field %s",field);

  if(flagId) {
    // Init for += syntax
    timeId = atom->timeId;
    if(SetFloatField(&timeId,arg,FieldPositive)==NULL) return(NULL);
  }
  else {