stft_tabulate.c

LastWave

    GGR = tab1GGR[i] = NewSignal();
    GGI = tab1GGI[i] = NewSignal();
    SizeSignal(GGR,GABOR_MAX_FREQID,YSIG);
    SizeSignal(GGI,GABOR_MAX_FREQID,YSIG);
    ZeroSig(GGR);
    ZeroSig(GGI);

    /* Fft */
    Fft1(shiftedWindow,NULL,GGR,GGI,1,NO);

    // The zero frequency : this must be strictly ensured
    GGR->Y[0] = realGG = 1.0;
    GGI->Y[0] = imagGG = 0.0;
    // The other frequencies
    for(freqId = 1; freqId<=GABOR_MAX_FREQID/4 && realGG*realGG+imagGG*imagGG>precision; freqId++) {
      realGG = GGR->Y[2*freqId];
      imagGG = GGI->Y[2*freqId];
      // Deal with computation round-off errors : 
      // -if the window has symmetry the imagGG should be zero
      // -in any case the resulting (realGG,imagGG) should have modulus at most 1
      if (flagAsymWindow==NO) imagGG = 0.0;
      if(realGG>1.0) {
	realGG=1.0;
	imagGG=0.0;
      }
      if(realGG<-1.0) {
	realGG=-1.0;
	imagGG=0.0;
      }
      // This situation has only been observed so far with spline3,exponential and FoF windows
      // for windowSize=16,32,64 and freqId = 1,2,3 with excess of modulus of the order of 1e-07
      if(realGG*realGG+imagGG*imagGG>1.0) {
	if(realGG*realGG+imagGG*imagGG-1.0>1e-06)
	  Errorf("StftTabulateGG : (WEIRED) modulus of GG exceeds 1 by more than %g",1e-06);
	imagGG=0.0;
      }
      GGR->Y[freqId] = realGG;
      GGI->Y[freqId] = imagGG;
    }
    for(;freqId<GGR->size; freqId++) {
      GGR->Y[freqId] = 0.0;
      GGI->Y[freqId] = 0.0;
    }
  }
  // DEBUG
  Printf("\n");
  DeleteSignal(paddedWindow);
  DeleteSignal(shiftedWindow);
}

// Gets a GG from the the tabulation
// -If the tabulation was not done for the windowShape, then we do it for all windowSizes
// -If the windowSize is not tabulated (i.e. so far not a power of 2), an ERROR is generated
// -TODO should we COMPUTE the GGR ? But should we allocate it ?????
void GetTabGG(char windowShape, unsigned long windowSize,SIGNAL *GGRptr,SIGNAL *GGIptr)
{
  SIGNAL *tab1GGR  = NULL;
  SIGNAL *tab1GGI  = NULL;
  unsigned short i;
  /* Check args */
  if(!WindowShapeIsOK(windowShape)) Errorf("GetTabGG : bad windowShape %d",windowShape);
  if(GGRptr==NULL||GGIptr==NULL)    Errorf("GetTabGG : NULL output");

  *GGRptr = NULL;
  *GGIptr = NULL;
  // If the GGR+GGI were not tabulated we do it now
  if(stftTabGGR[windowShape]==NULL) StftTabulateGG(windowShape);
  tab1GGR = stftTabGGR[windowShape];
  tab1GGI = stftTabGGI[windowShape];
  if(tab1GGR==NULL) Errorf("GetTabGG: NULL tab1GGR");
  if(tab1GGI==NULL) Errorf("GetTabGG: NULL tab1GGI");


  // Find the content at the desired windowSize
  i = GetIndexWindowSize(windowSize);
  if(i<nTabWindowSizes) {
    *GGRptr = tab1GGR[i];
    *GGIptr = tab1GGI[i];
  } else   
    Errorf("GetTabGG : windowSize = %d is not tabulated",windowSize);
  // DEBUG
  if(*GGRptr==NULL) Errorf("GetTabGG: NULL realGG");
  if(*GGIptr==NULL) Errorf("GetTabGG: NULL imagGG");
}


// stftTabWindowShapes and stftTabWindowFactors are initialised in InitStftTabulations
static SIGNAL** stftTabWindowShapes  = NULL;
static LWFLOAT**  stftTabWindowFactors = NULL;
// If the windows with the given 'windowShape' are not tabulated,
// we load them from a file or tabulate them,
// else just do nothing.
static void StftTabulateWindowShape(char windowShape)
{
  LWFLOAT (*f)(SIGNAL,LWFLOAT);

  unsigned short i;
  SIGNAL* windowShapes = NULL;
  LWFLOAT*  windowFactors= NULL;

  /* Checking arguments */
  if(!WindowShapeIsOK(windowShape)) 
    Errorf("StftTabulateWindowShape : bad windowShape %d",windowShape);
  // Quit without doing anything if already tabulated
  if(stftTabWindowShapes[windowShape]!=NULL) return;

  // TODO ?? : try to load from file here

  // Getting the right window
  GetWindowShapeFunc(windowShape,&f);

  /* Allocating */
  stftTabWindowShapes[windowShape]  = (SIGNAL*) Calloc(nTabWindowSizes,sizeof(SIGNAL));
  stftTabWindowFactors[windowShape] = (LWFLOAT *) Calloc(nTabWindowSizes,sizeof(LWFLOAT));
  windowShapes  = stftTabWindowShapes[windowShape];
  windowFactors = stftTabWindowFactors[windowShape];

  // Allocating the signals and filling them
  for(i=0; i<nTabWindowSizes; i++) {
    windowShapes[i] = NewSignal();
    SizeSignal(windowShapes[i],stftTabWindowSizes[i],YSIG);
    windowFactors[i]=(*f)(windowShapes[i],0.0);
  }
}

// Gets a window from the the tabulated windows.
// -If the tabulation was not done for the windowShape, then we do it for all windowSizes
// -If the windowSize is not tabulated (i.e not a power of 2) an ERROR is generated
// - TODO should we COMPUTE the window ? But should we allocate it ?????
void GetTabWindow(char windowShape,unsigned long windowSize,SIGNAL *windowPtr)
{
  unsigned short i;

  /* Checking arguments */
  if(!WindowShapeIsOK(windowShape)) 
    Errorf("GetTabWindow : bad windowShape %d",windowShape);
  if(windowPtr == NULL)
    Errorf("GetTabWindow : NULL output");
    
  // If the windowShape was not tabulated we do it now
  if(stftTabWindowShapes[windowShape]==NULL) 
    StftTabulateWindowShape(windowShape);

  // Find the tabulated windowShape at the right 'windowSize'
  i=GetIndexWindowSize(windowSize);
  if(i<nTabWindowSizes) {
    *windowPtr = stftTabWindowShapes[windowShape][i];
    // DEBUG
    if((*windowPtr)->size!=windowSize) 
      Errorf("GetTabWindow : corrupted windows");
    return;
  } else
    Errorf("GetTabWindow : windowSize %d is not tabulated",windowSize);
}

// Gets a windowFactor from the the tabulated windows.
// -If the tabulation was not done for the windowShape, then we do it for all windowSizes
// -If the windowSize is not tabulated (i.e not a power of 2) an ERROR is generated
// - TODO should we COMPUTE the window ? But should we allocate it ?????
void GetTabWindowFactor(char windowShape,unsigned long windowSize,LWFLOAT *windowFactorPtr)
{
  unsigned short i;

  /* Checking arguments */
  if(!WindowShapeIsOK(windowShape)) 
    Errorf("GetTabWindowFactor : bad windowShape %d",windowShape);
  if(windowFactorPtr == NULL)
    Errorf("GetTabWindowFactor : NULL output");
    
  // If the windowShape was not tabulated we do it now
  if(stftTabWindowShapes[windowShape]==NULL) 
    StftTabulateWindowShape(windowShape);

  // Find the tabulated windowShape at the right 'windowSize'
  i=GetIndexWindowSize(windowSize);
  if(i<nTabWindowSizes) {
    *windowFactorPtr = stftTabWindowFactors[windowShape][i];
    return;
  } else
    Errorf("GetTabWindowFactor : windowSize %d is not tabulated",windowSize);
}

//
// The signals where complex exponentials are tabulated
//
static SIGNAL expikR=NULL;
static SIGNAL expikI=NULL;

//  Function to tabulate complex exponentials
//  TODO ?? - try to load them from a file
//  Else we allocate and tabulate them
static void StftTabulateComplexExponentials(void)
{
  unsigned long i;
  // TODO ?? : try to load from file here

  //
  // Allocating and tabulating the finest resolution complex exponentials
  //
  expikR = NewSignal();
  expikI = NewSignal();
  SizeSignal(expikR,GABOR_MAX_FREQID,YSIG);
  SizeSignal(expikI,GABOR_MAX_FREQID,YSIG);
  ZeroSig(expikR);
  ZeroSig(expikI);

  // Tabulation of the frequency 0
  expikR->Y[0] = 1.0;
  expikI->Y[0] = 0.0;
  // Tabulation of the Nyquist frequency
  expikR->Y[GABOR_NYQUIST_FREQID] = -1.0;
  expikI->Y[GABOR_NYQUIST_FREQID] = 0.0;
  // Tabulation of the other frequencies
  for(i=1; i < GABOR_NYQUIST_FREQID; i++) {
    expikR->Y[GABOR_MAX_FREQID-i] = expikR->Y[i] = cos(2.0*(M_PI*i)/GABOR_MAX_FREQID);
    expikI->Y[i]            = sin(2.0*(M_PI*i)/GABOR_MAX_FREQID);
    expikI->Y[GABOR_MAX_FREQID-i] = -expikI->Y[i];
  }
}

// Gets the complex exponentials from the tabulation.
void GetTabExponentials(SIGNAL *expikRptr,SIGNAL *expikIptr)
{
  /* Checking arguments */
  if(expikRptr==NULL||expikIptr==NULL) Errorf("GetTabExponentials : NULL output");
  *expikRptr = expikR;
  *expikIptr = expikI;
}

// This function is called when the Stft package is loaded
void InitStftTabulations(void)
{
  unsigned long windowSize;
  unsigned long i;

  // Compute the number of values of 'windowSize' for which tabulation occurs
  for(nTabWindowSizes=0,windowSize=4; windowSize <= STFT_MAX_WINDOWSIZE; nTabWindowSizes++, windowSize<<=1) {}
  // Allocating an computing the map i -> windowSize
  stftTabWindowSizes = (unsigned long*) Calloc(nTabWindowSizes,sizeof(unsigned long));
  for(i=0,windowSize=4; i<nTabWindowSizes; i++, windowSize<<=1) {
    stftTabWindowSizes[i]=windowSize;
  }

  // Allocate and tabulate all complex exponentials
  StftTabulateComplexExponentials();

  // Allocate memory for future tabulation of windows for every windowShape,
  // Allocation for a given windowShape will only occur if necessary.
  stftTabWindowShapes  = (SIGNAL **)  Calloc(LastWindowShape,sizeof(SIGNAL*));
  stftTabWindowFactors = (LWFLOAT **)   Calloc(LastWindowShape,sizeof(LWFLOAT*));
  stftTabGGR           = (SIGNAL **) Calloc(LastWindowShape,sizeof(SIGNAL*));
  stftTabGGI           = (SIGNAL **) Calloc(LastWindowShape,sizeof(SIGNAL*));
}

// This function is called when the Stft package is unloaded
// (for example when LastWave quits)
void FreeStftTabulations(void)
{
  unsigned long i;
  char windowShape;
  // Free memory allocated for tabulation of windows for every windowshape
  if(stftTabWindowShapes != NULL) {
    for(windowShape = 0; windowShape < LastWindowShape; windowShape++) {
      if(stftTabWindowShapes[windowShape] != NULL) {
	for(i=0; i < nTabWindowSizes; i++) {
	  if(stftTabWindowShapes[windowShape][i]!=NULL) {
	    DeleteSignal(stftTabWindowShapes[windowShape][i]);
	    stftTabWindowShapes[windowShape][i]=NULL;
	  }
	  else 
	    Warningf("FreeStftTabulations - Weird error - no signal found");
	}
	Free(stftTabWindowShapes[windowShape]);
	stftTabWindowShapes[windowShape] = NULL;
      }
    }
    Free(stftTabWindowShapes);
    stftTabWindowShapes = NULL;
  }
  else
    Warningf("FreeStftTabulations - stftTabWindowShapes is NULL!");
  // Free memory allocated for tabulation of real GG 
  if(stftTabGGR != NULL) {
    for(windowShape = 0; windowShape < LastWindowShape; windowShape++) {
      if(stftTabGGR[windowShape] != NULL) {
	for(i = 0; i < nTabWindowSizes; i++)
	  if(stftTabGGR[windowShape][i] != NULL) {
	    DeleteSignal(stftTabGGR[windowShape][i]);
	    stftTabGGR[windowShape][i] = NULL;
	  }
	Free(stftTabGGR[windowShape]);
	stftTabGGR[windowShape] = NULL;
      }
    }
    Free(stftTabGGR);
    stftTabGGR = NULL;
  }
  else
    Errorf("FreeStftTabulations - stftTabGGR is NULL!");

  //  Free memory allocated for tabulation of imag GG
  if(stftTabGGI != NULL) {
    for(windowShape = 0; windowShape < LastWindowShape; windowShape++) {
      if(stftTabGGI[windowShape] != NULL) {
	for(i = 0; i < nTabWindowSizes; i++)
	  if(stftTabGGI[windowShape][i] != NULL) {
	    DeleteSignal(stftTabGGI[windowShape][i]);
	    stftTabGGI[windowShape][i] = NULL;
	  }
	Free(stftTabGGI[windowShape]);
	stftTabGGI[windowShape] = NULL;
      }
    }
    Free(stftTabGGI);
    stftTabGGI = NULL;
  }
  else
    Warningf("FreeStftTabulations - stftTabGGI is NULL!");

  Free(stftTabWindowSizes);
  stftTabWindowSizes = NULL;
  nTabWindowSizes    = 0;
}
/* EOF */