stft_window.c

LastWave

  int j;
  LWFLOAT limit;
  LWFLOAT i;
  LWFLOAT beta;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  LWFLOAT a;
  LWFLOAT expon;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("FoF : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("FoF : bad shift %g",shift);
  
  /* Damping factor */
  
  a = log(decayFoF);
  
  
  /* scale */
  
  expon = a/window->size; 
  beta = betaFoF/window->size;
  
  /* limit of the cos */
  
  limit = M_PI/beta;
  
  /* The window */
  
  for(j=0, i=0.0; i <= limit ; j++, i += 1.0) {
    window->Y[j] = 0.5*(1-cos(beta*i))*exp(-expon*i);
    energy += window->Y[j]*window->Y[j];
  }
  
  for(; i < window->size; j++, i += 1.0) {
    window->Y[j] = exp(-expon*i);
    energy += window->Y[j]*window->Y[j];
  }
  
  /* Putting the first point to zero*/
  
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  
  factor = 1/sqrt(energy);
  
  for(j = 0; j < window->size; j++) {
    window->Y[j] *= factor;
  }
  /*  Printf("FoF:(a: %f,factor: %f,size: %d,limit: %f)\n",a,factor,window->size,limit);  javi */
  
  return(factor);
}

/*************************/
/*    Asymmetric3        */
/*************************/

/*  Another asymmetric window */


static LWFLOAT Asymmetric3(SIGNAL window,LWFLOAT shift)
{
  LWFLOAT decay = 1e4;  
  int j;
  LWFLOAT limit;
  LWFLOAT beta = M_PI/(log(window->size)/log(2)-1);	/* beta = PI/octave */
  LWFLOAT i;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  LWFLOAT a;
  LWFLOAT expon;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Asymmetric : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Asymmetric : bad shift %g",shift);
  
  /* Damping factor */
  
  a = log(decay);
  expon = a/window->size;  /* scaled */
  
  /* limit of the cos */
  
  limit = M_PI/beta;
  
  /* The window */
  
  for(j=0, i=0.0; i <= limit ; j++, i += 1.0) {
    window->Y[j] = 0.5*(1-cos(beta*i))*exp(-expon*limit);
    energy += window->Y[j]*window->Y[j];
  }
  
  for(; i < window->size; j++, i += 1.0) {
    window->Y[j] = exp(-expon*i);
    energy += window->Y[j]*window->Y[j];
  }
  
  /* Putting the first point to zero*/
  
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  
  factor = 1/sqrt(energy);
  
  for(j = 0; j < window->size; j++) {
    window->Y[j] *= factor;
  }
  /*  Printf("Asymmetric Window:(a: %f,factor: %f,size: %d,limit: %f)\n",a,factor,window->size,limit);  javi */
  
  return(factor);
}

/*  asymmetric window */
static LWFLOAT Asymmetric32(SIGNAL window,LWFLOAT shift)
{
  LWFLOAT decay = 1e4;  
  int j;
  LWFLOAT i;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  LWFLOAT a;
  LWFLOAT expon;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Asymmetric : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Asymmetric : bad shift %g",shift);
  
  /* Damping factor */
  
  a = log(decay);
  expon = a/window->size;  /* scaled */
  
  
  /* The window */
  
  for(j=0, i=0; i < window->size; j++, i++) {
    window->Y[j] = (i+1)*exp(-expon*i);
    energy += window->Y[j]*window->Y[j];
  }
  
  
  /* Putting the first point to zero*/
  
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  
  factor = 1/sqrt(energy);
  
  for(j = 0; j < window->size; j++) {
    window->Y[j] *= factor;
  }
  Printf("Asymmetric Window:(a: %f,factor: %f,size: %d)\n",a,factor,window->size); /* javi */
  
  return(factor);
}


/**********************************************/
/*      Get the window function            */
/**********************************************/

void GetWindowShapeFunc(char windowShape, LWFLOAT (**f)(SIGNAL,LWFLOAT))
{
  if(!WindowShapeIsOK(windowShape))
    Errorf("GetWindowShapeFunc : unknown window shape %d",windowShape);
  
  switch(windowShape) {
  case Spline0WindowShape :
    *f = &Spline0;
    break;
  case Spline1WindowShape :
    *f = &Spline1; 
    break;
  case Spline2WindowShape : 
    *f = &Spline2; 
    break;
  case Spline3WindowShape : 
    *f = &Spline3; 
    break;
  case GaussWindowShape : 
    *f = &Gauss; 
    break;
  case HammingWindowShape : 
    *f = &Hamming; 
    break;
  case HanningWindowShape : 
    *f = &Hanning; 
    break;
  case BlackmanWindowShape : 
    *f = &Blackman; 
    break;
  case ExponentialWindowShape : 
    *f = &Exponential; 
    break;
  case FoFWindowShape : 
    *f = &FoF; 
    break;
  case Asym3WindowShape : 
    *f = &Asymmetric3; 
    break;
  default: 
    Errorf("GetWindowShapeFunc : bad window shape %d",windowShape);
    break;
  }
}

char  Name2WindowShape(char *name)
{
  if(!strcmp(name,"spline0")) return(Spline0WindowShape);
  if(!strcmp(name,"spline1")) return(Spline1WindowShape);
  if(!strcmp(name,"spline2")) return(Spline2WindowShape);
  if(!strcmp(name,"spline3")) return(Spline3WindowShape);
  if(!strcmp(name,"gauss")) return(GaussWindowShape);
  if(!strcmp(name,"hamming")) return(HammingWindowShape);
  if(!strcmp(name,"hanning")) return(HanningWindowShape);
  if(!strcmp(name,"blackman")) return(BlackmanWindowShape);
  if(!strcmp(name,"exponential")) return(ExponentialWindowShape);
  if(!strcmp(name,"FoF")) return(FoFWindowShape);
  if(!strcmp(name,"asymmetric3")) return(Asym3WindowShape);
  Errorf("Unknown window name '%s'",name);
  
  // This should never be reached but the compiler may want to be sure the function returns a value
  return(LastWindowShape);
}

char * WindowShape2Name(char windowShape)
{
  if(!WindowShapeIsOK(windowShape))
    Errorf("WindowShape2Name : Unknown window shape %d",windowShape);
  
  switch(windowShape) {
  case Spline0WindowShape : return("spline0");
  case Spline1WindowShape : return("spline1");
  case Spline2WindowShape : return("spline2");
  case Spline3WindowShape : return("spline3");
  case GaussWindowShape   : return("gauss");
  case HammingWindowShape   : return("hamming");
  case HanningWindowShape   : return("hanning");
  case BlackmanWindowShape   : return("blackman");
  case ExponentialWindowShape   : return("exponential");
  case FoFWindowShape   : return("FoF");
  case Asym3WindowShape   : return("asymmetric3");	
  }
  Errorf("Unknown window shape %d",windowShape);
  // This should never be reached but the compiler may want to be sure the function returns a value
  return("unknown");
}

#ifdef STFT_ADVANCED
void C_GetWindowShape(char **argv)
{
  char *windowShapeName;
  LWFLOAT (*f) (SIGNAL,LWFLOAT);
  
  SIGNAL signal;
  LWFLOAT shift;
  int size;
  LWFLOAT factor;
  
  argv = ParseArgv(argv,tSIGNAL,&signal,tINT,&size,tFLOAT_,0.0,&shift,tSTR_,"gauss",&windowShapeName,0);
  if(shift < 0 || shift >= 1)
    Errorf("shift = %g should remain within [0,1[");
  GetWindowShapeFunc(Name2WindowShape(windowShapeName),&f);
  SizeSignal(signal,size,YSIG);
  factor = (*f)(signal,shift);
  /*    SetResultFloat(factor); */
  return;
}
#endif

char GetFlagAsymWindowShape(char windowShape)
{
  char flagAsym;
  
  switch(windowShape) {
  case ExponentialWindowShape :
    flagAsym = YES;
    break;
  case FoFWindowShape :
    flagAsym = YES;
    break;
  case Asym3WindowShape :
    flagAsym = YES;
    break;
  default:
    flagAsym = NO;
    break;
  }
  return(flagAsym);
}

/* 
 * 'Max window' is the index of the 'maximum'
 * of the window (sometimes it is not actually the maximum ...
 * in the array of size 'size', 
 * where indexes run from 0 to size-1
 */
int GetMaxWindowShape(char windowShape,int size)		/* javi */
{
  int max;
  
  switch(windowShape) {	
  case ExponentialWindowShape:
    max = 1;
    break;
  case FoFWindowShape:
    /*	max = GetMaxFoF(size); */
    max = 1;
    break;
  case Asym3WindowShape:
    max = (int) (log(size)/log(2)-1);
    break;
  default:
    max = size/2;
    break;
  }
  return(max);
}

int GetMaxFoF(int size)
{
  int Max;
  int octave;
  
  octave = (int) (log(size)/log(2));
  
  switch(octave) {	
  case 2:
    Max = 1;
    break;
  case 3:
    Max = 1;
    break;
  case 4:
    Max = 1;
    break;
  case 5:
    Max = 3;
    break;
  case 6:
    Max = 5;    
    break;
  case 7:
    Max = 10;
    break;
  case 8:
    Max = 21;
    break;
  case 9:
    Max = 42;
    break;
  case 10:
    Max = 84;
    break;
  case 11:
    Max = 168;
    break;
  case 12:
    Max = 335;
    break;
  case 13:
    Max = 670;
    break;
  case 14:
    Max = 1341;
    break;
  case 15:
    Max = 2682;
    break;
  case 16:
    Max = 5364;
    break;
  default:
    Errorf("GetMaxFof : octave %d not treated",octave);
    Max = 1; // The compiler wants to see an initialized return value
    break;
  }
  return(Max);
}

/* Computes the time support of a window, that is to say 
 * the closed interval [timeIdMin,timeIdMax] on which the
 * window is non zero
 * WARNING : it is the calling function's responsibility
 *           to check whether this support intersects the
 *           range of indexes of an array (e.g. [0,signalSize-1])
 */
void ComputeWindowSupport(int windowSize,char windowShape,LWFLOAT timeId,
			  int *pTimeIdMin,int *pTimeIdMax)
{
  int maxWindowShape;
  
  /* Checking arguments */
  if(!INRANGE(1,windowSize,STFT_MAX_WINDOWSIZE))
    Errorf("ComputeWindowSupport : bad windowSize %d",windowSize);
  if(!WindowShapeIsOK(windowShape))
    Errorf("ComputeWindowSupport : unknown windowShape %d",windowShape);
  if(pTimeIdMin == NULL || pTimeIdMax == NULL)
    Errorf("ComputeWindowSupport : NULL output");
  
  /* Computing */
  // Location index of the maximum of the window
  maxWindowShape=GetMaxWindowShape(windowShape,windowSize);
  
  if(timeId == (int) timeId) {
    /* Checked correct on the 20/02/2001 by R. Gribonval
     *
     * a/timeId-maxWindowShape is the index of the first 
     * point of the window, which value is set to zero, hence the +1
     * b/the indexes range from 0 to windowSize-1, shifted by 
     * timeId-maxWindowShape
     */
    *pTimeIdMin = (int) timeId-maxWindowShape+1;
    *pTimeIdMax = (int) timeId-maxWindowShape+windowSize-1; 
  }
  /* We should define what happens when timeId != integer */
  else
    Errorf("ComputeWindowSupport : timeId %g is not an integer, behavior undefined",timeId);
  
  /* *pTimeIdMin = timeId-windowSize2+1;
   *pTimeIdMax = timeId+windowSize2-1;  */
}

/* EOF */