stft.c

LastWave

  if(stft->fRate*(stft->nSubBands-1) != GABOR_NYQUIST_FREQID)
    Errorf("CheckStft %d %s : bad freq grid %d * %d != %d = Nyquist",stft->windowSize,StftType2Name(stft->type),stft->fRate,stft->nSubBands-1,GABOR_NYQUIST_FREQID);
  /* Type dependent checkings */
  if(stft->signalSize != 0) {
    switch(stft->type) {
    case ComplexStft:
      if(stft->real == NULL || stft->imag == NULL)
	Errorf("CheckStft %d %s : NULL complex data",stft->windowSize,StftType2Name(stft->type));
      break;
    case RealStft:
    case PhaseStft:
    case HighResStft:
    case HarmoStft:
      if(stft->real == NULL)
	Errorf("CheckStft %d %s : NULL real data",stft->windowSize,StftType2Name(stft->type));
      if(stft->imag != NULL)
	Errorf("CheckStft %d %s : NON NULL imag data",stft->windowSize,StftType2Name(stft->type));
      break;
    default :
      Errorf("CheckStft : (Weird) unknown type %d",stft->type);
    }
  }
}

void CheckStftNotEmpty(const STFT stft)
{
  CheckStft(stft);
  if(stft->signalSize == 0)  Errorf("CheckStftNotEmpty : empty stft");
}

void CheckStftGridCompat(const STFT stftGrid,const STFT stftSubGrid)
{
  CheckTFContentCompat(stftGrid,stftSubGrid);
  if(stftGrid->windowSize     != stftSubGrid->windowSize ||
     stftGrid->borderType != stftSubGrid->borderType ||
     stftGrid->windowShape != stftSubGrid->windowShape)
    Errorf("CheckStftGridCompat : Weird error");
  if(stftGrid->tRate % stftSubGrid->tRate != 0)
    Errorf("CheckStftGridCompat : incompatible tRates %d %d",stftGrid->tRate,stftSubGrid->tRate);
  if(stftGrid->fRate % stftSubGrid->fRate != 0)
    Errorf("CheckStftGridCompat : incompatible fRates %d %d",stftGrid->fRate,stftSubGrid->fRate);
}

void CheckStftComplex(const STFT stft)
{
  CheckStftNotEmpty(stft);
  if(stft->type != ComplexStft)
    Errorf("CheckStftComplex : stft is not complex");
}

void CheckStftReal(const STFT stft)
{
  CheckStftNotEmpty(stft);
  if(stft->type != RealStft)
    Errorf("CheckStftReal : stft is not energy");
}

void CheckStftPhase(const STFT stft)
{
  CheckStftNotEmpty(stft);
  if(stft->type != PhaseStft)
    Errorf("CheckStftPhase : stft is not energy");
}

void CheckStftHighRes(const STFT stft)
{
  CheckStftNotEmpty(stft);
  if(stft->type != HighResStft)
    Errorf("CheckStftHighRes : stft is not highres");
}

void CheckStftHarmo(const STFT stft)
{
  CheckStftNotEmpty(stft);
  if(stft->type != HarmoStft)
    Errorf("CheckStftHarmo : stft is not harmo");
}

void CheckStftUpToDate(const STFT stft)
{
  CheckStftNotEmpty(stft);
  if(stft->flagUpToDate == NO)
    Errorf("CheckStftUpToDate : stft is out of date");
}


/* Allocation/Desallocation of STFT */
static void InitStft(STFT stft)
{	
  InitTFContent(stft);

  stft->tRate      = 1;
  stft->nFrames    = 0;
  stft->fRate      = 1;
  stft->nSubBands    = 0;
  stft->borderType = STFT_DEFAULT_BORDERTYPE;
  stft->firstp     = 0;
  stft->lastp      = -1;
  stft->freqIdMin  = 0;
  stft->freqIdMax  = 0;
  stft->windowShape = STFT_DEFAULT_WINDOWSHAPE;
  stft->windowSize 	   = 0;
  stft->flagUpToDate 	= NO;
  
  /* Default type for a stft : complex, zero-padded, gabor window */
  stft->type		= STFT_DEFAULT_STFTTYPE;
  /* The complex part */
  stft->real 		= NULL;
  stft->imag 		= NULL;

  /* The harmo part */
  stft->iFMO = 1;
  stft->dimHarmo             = 1;
}

STFT NewStft(void)
{
  STFT stft;
  
#ifdef DEBUGALLOC
  DebugType = "Stft";
#endif
  
  stft = (STFT) Malloc(sizeof(Stft));
  InitValue(stft,&tsStft);
  InitStft(stft);
  return(stft);
}

void ClearStft(STFT stft)
{
  if(stft == NULL)
    Errorf("ClearStft : NULL stft");
  
  if(stft->real) {
    DeleteSignal(stft->real);
    stft->real = NULL;
  }
  if(stft->imag) {
    DeleteSignal(stft->imag);
    stft->imag = NULL;
  }

  InitStft(stft);
}

STFT DeleteStft(STFT stft)
{
  if(stft == NULL)
    Errorf("DeleteStft : NULL stft");
  
  if (stft->nRef==0) {
    Errorf("*** Danger : trying to delete a temporary stft\n");
  }
  
  RemoveRefValue(stft);
  if (stft->nRef > 0) 
    return(NULL);
  
  if(stft->real) {
    DeleteSignal(stft->real);
    stft->real = NULL;
  }
  if(stft->imag) {
    DeleteSignal(stft->imag);
    stft->imag = NULL;
  }
  
#ifdef DEBUGALLOC
  DebugType = "Stft";
#endif
  
  Free(stft);
  return(NULL);
}

void TouchStft(STFT stft)
{
  if(stft == NULL)
    Errorf("TouchStft : NULL stft");
  stft->flagUpToDate = NO;
}



void PrintStft(const STFT stft,char flagShort)
{
  /* Checking argument */
  CheckStft(stft);
  
  if(flagShort) {
    Printf("[%3d]    %-7d %-7d %-7d %-7d %s %s\n",
	   stft->windowSize,
	   stft->tRate,stft->nFrames,stft->fRate,stft->nSubBands,
	   BorderType2Name(stft->borderType),WindowShape2Name(stft->windowShape),StftType2Name(stft->type));
  }
  else {
    /* Case of an empty stft */
    if(stft->windowSize <= 1) {
      Printf("  empty\n");
      return;
    }
    Printf("  window        : %s %d\n",
	   WindowShape2Name(stft->windowShape),stft->windowSize);
    Printf("  type          : %s\n",StftType2Name(stft->type));
    Printf("  size          : [%d frames x %d subbands]\n",stft->nFrames,stft->nSubBands);
    Printf("  tRate,fRate   : %-4d %-4d\n",stft->tRate,stft->fRate);
    if(stft->type==HarmoStft) {
      Printf("  range of f0   : [%g %g] ([%d %d])\n",
	     FreqId2Freq(stft,stft->freqIdMin),FreqId2Freq(stft,stft->freqIdMax),stft->freqIdMin,stft->freqIdMax);
      Printf("  dimension     : %d\n",stft->dimHarmo);
      Printf("  quasi-ortho   : %d\n",stft->iFMO);
    }
    Printf("Analyzed signal information\n");
    Printf("  border     : %s\n",BorderType2Name(stft->borderType));
    Printf("  signalSize : %d\n",stft->signalSize);
    Printf("  x0         : %g\n",stft->x0);
    Printf("  dx         : %g [Sampling Frequency = %g]\n",stft->dx,1/stft->dx);
    Printf("  firstp     : %d (%g)\n",stft->firstp,TimeId2Time(stft,stft->firstp));
    Printf("  lastp      : %d (%g)\n",stft->lastp,TimeId2Time(stft,stft->lastp));
  }
}

/*********************************************************/
/*
 *	Various functions to set a stft field
 */
/*********************************************************/





/*
 * Compute the frequency 'box' [freqIdMin,freqIdMax] on which 
 * one can look for partial 0<=k<stft->dimHarmo corresponding to 
 * the fundamental frequency freq0Id.
 */
char HarmoPartialBox(STFT stft,int k,int freq0Id, 
		     int* pk_f0IdMin,int* pk_f0IdMax)
{  
  LWFLOAT A = 0;
  LWFLOAT df;
  /* Checking arguments */
  CheckStftHarmo(stft); 
  if(k<0 || k >= stft->dimHarmo)
    Errorf("HarmoPartialBox : k %d is out of range [0, %d[",k,stft->dimHarmo);

  *pk_f0IdMax = (k+1)*freq0Id;
  *pk_f0IdMin = (k+1)*freq0Id;
  
  /* Not really used for now since A = 0 */
  df = (int) ((A/(stft->windowSize*stft->dx))*stft->dx*GABOR_MAX_FREQID+.5);
  
  *pk_f0IdMin -= (int) df; 
  *pk_f0IdMax += (int) df;
  
  if(*pk_f0IdMin > *pk_f0IdMax) {
    Errorf("k_f0IdMin > k_f0IdMax\n");
  }
  
  if(*pk_f0IdMin >= stft->fRate*stft->nSubBands) {
    //#define DEBUGPARTIAL
#ifdef DEBUGPARTIAL
    Printf("Bad Box (1): s=%d k=%d f0=%d",stft->windowSize,k,freq0Id);
    Printf("box=[%d %d] rate=%d length=%d\n",*pk_f0IdMin,*pk_f0IdMax,stft->fRate,stft->nSubBands);
#endif
    return(NO);
  }
  
  if(*pk_f0IdMax > stft->fRate*(stft->nSubBands-1)) {
#ifdef DEBUGPARTIAL
    Printf("Bad Box (2): s=%d k=%d f0=%d",stft->windowSize,k,freq0Id);
    Printf("box=[%d %d] rate=%d length=%d\n",*pk_f0IdMin,*pk_f0IdMax,stft->fRate,stft->nSubBands);
#endif
    return(NO);
  }
#ifdef DEBUGPARTIAL
  if(k != 0) {
    Printf("Good Box (2): o=%d k=%d f0=%d",stft->windowSize,k,freq0Id);
    Printf("box=[%d %d] rate=%d length=%d\n",*pk_f0IdMin,*pk_f0IdMax,stft->fRate,stft->nSubBands);
  }
#endif
  return(YES);
  
  /* Second version NEVER REACHED, NOT FINALIZED */
  *pk_f0IdMin = (k+1)*freq0Id-stft->fRate*((freq0Id/stft->fRate - stft->iFMO)/2);
  
  if(*pk_f0IdMin >= stft->fRate*stft->nSubBands) 
    return(NO);
  
  *pk_f0IdMax = (k+1)*freq0Id+stft->fRate*((freq0Id/stft->fRate - stft->iFMO)/2);
  *pk_f0IdMax = MIN(*pk_f0IdMax-stft->fRate,stft->fRate*(stft->nSubBands-1));
  
  if(*pk_f0IdMin > *pk_f0IdMax) 
    Errorf("k_f0IdMin > k_f0IdMax\n");
  
  return(YES);
  
}

/**********************************************************************/
/*
 * Compute a time-freq grid at a given windowSize, given two numbers
 * 'timeRedund' and 'freqRedund' that specify whether there is subsampling
 * or not at the given scale.
 */
/**********************************************************************/
void SizeStftGrid(STFT stft,int windowSize,int timeRedund,int freqRedund)
{
  CheckTFContent(stft);
  if(stft->signalSize == 0)  Errorf("SizeStftGrid : bad signalSize %d <= 0",stft->signalSize);
  /* Is it necessary ? */
  if(!INRANGE(STFT_MIN_WINDOWSIZE,windowSize,STFT_MAX_WINDOWSIZE) || (windowSize > GABOR_MAX_FREQID))
    Errorf("SizeStftGrid : Bad windowSize %d not in [%d %d] or larger that %d",windowSize,STFT_MIN_WINDOWSIZE,STFT_MAX_WINDOWSIZE,GABOR_MAX_FREQID);
  // TODO : understand or remove ?
  if(!INRANGE(1,timeRedund,GABOR_MAX_FREQID))
    Errorf("SizeStftGrid : Bad timeRedund %d for GABOR_MAX_FREQID = %d",timeRedund,GABOR_MAX_FREQID);
  if(!INRANGE(1,freqRedund,GABOR_MAX_FREQID))
    Errorf("SizeStftGrid : Bad freqRedund %d for GABOR_MAX_FREQID = %d",freqRedund,GABOR_MAX_FREQID);

  // In the time domain ... use of timeRedund
  if (windowSize <= timeRedund)    stft->tRate 	= 1;
  else                             stft->tRate 	= windowSize/(timeRedund);

  // In the frequency domain ... use of freqRedund !
  if (windowSize > GABOR_MAX_FREQID/freqRedund)
    stft->fRate 	= 1;
  else
    stft->fRate 	= 2*GABOR_MAX_FREQID/((freqRedund)*windowSize);

  stft->nFrames	= stft->signalSize/(stft->tRate);
  if((stft->nFrames)*(stft->tRate)< stft->signalSize)
    (stft->nFrames) += 1;

  stft->nSubBands 	= GABOR_NYQUIST_FREQID/(stft->fRate)+1;
  stft->windowSize = windowSize;
}

/***********************************************/
/*
 * Setting a stft time-frequency grid
 * and allocating the data arrays.
 * Keeps the ATFCONTENT fields
 * Resets it to STFT_DEFAULT_BORDERTYPE and STFT_DEFAULT_WINDOWSHAPE
 */
/***********************************************/
void SizeStft(STFT stft,char windowShape,int windowSize,
	      int timeRedund,int freqRedund,
	      char borderType,char stftType,
	      LWFLOAT freq0IdMin,LWFLOAT freq0IdMax,
	      int iFMO,int dimHarmo)
{
  struct aTFContent tfContent;
  // Used for the grid
  int shiftMin,shiftMax;
  // Only for HarmoStft type
  int f0IdMin, f0IdMax;
  int i_f0Min,i_f0Max;
  int size;
  int octave;
  
  // Checking arguments
  if(!WindowShapeIsOK(windowShape))   Errorf("SizeStft : Bad windowShape %d\n",windowShape);
  // TODO REMOVE THIS CONSTRAINT !!!!
  octave = (int) (log(windowSize)/log(2.0)+0.5);
  if((1<<octave) != windowSize)       Errorf("SizeStft : windowSize %d is not a power of two!",windowSize);

  if(!BorderTypeIsOK(borderType))     Errorf("SizeStft : unknown border type %d",borderType);
  // Checkings specific to HarmoStft type
  if(stftType == HarmoStft) {
    if(freq0IdMin <= 0.0)
      Errorf("SizeStft : bad freq0IdMin %g <= 0.0",freq0IdMin);
    if(freq0IdMin > freq0IdMax)
      Errorf("SizeStft : bad freq0 range [%g %g]!",freq0IdMin,freq0IdMax);
    if(freq0IdMax > GABOR_NYQUIST_FREQID)
      Errorf("SizeStft : bad freq0IdMax %g > %d",freq0IdMax,GABOR_NYQUIST_FREQID);
    if(iFMO <= 0)
      Errorf("SizeStft : bad iFMO %d",iFMO);
    if(dimHarmo <= 0)
      Errorf("SizeStft : bad dimHarmo %d",dimHarmo);
  }
  
  
  // Clearing the signals if necessary
  CopyFieldsTFContent(stft,&tfContent);
  ClearStft(stft);
  CopyFieldsTFContent(&tfContent,stft);
  
  // Setting the time freq grid 
  SizeStftGrid(stft,windowSize,timeRedund,freqRedund);

  // The window and border now
  stft->windowShape = windowShape;
  stft->borderType = borderType;
  
  /* 
   * Determines which are the first and last timeId in a stft
   * that are not affected by left/right side effects
   */   
  
  /* 
   * The time support of a window centered at timeId is 
   * [timeId+shiftMin,timeId+shiftMax]
   * It is affected by border effects only if this support
   * intersects the outside of [0,signalSize-1], hence
   * lastp + shiftmax = signalSize-1
   * firstp+ shiftmin = 0
   */
  ComputeWindowSupport(stft->windowSize,stft->windowShape,0.0,
		       &shiftMin,&shiftMax);
  
  stft->firstp = MAX(0,-shiftMin);
  stft->lastp  = MIN(stft->signalSize-1,stft->signalSize-1 - shiftMax);
  TouchStft(stft);

  stft->freqIdMin = 0;
  stft->freqIdMax = GABOR_NYQUIST_FREQID;

  switch(stftType) {
  case ComplexStft:
    stft->type = stftType;
    if(stft->real == NULL)
      stft->real = NewSignal();
    if(stft->imag == NULL)
      stft->imag = NewSignal();
    SizeSignal(stft->real,stft->nFrames*stft->nSubBands,YSIG); 
    SizeSignal(stft->imag,stft->nFrames*stft->nSubBands,YSIG); 
    ZeroSig(stft->real);
    ZeroSig(stft->imag);
    stft->real->dx = stft->fRate;
    stft->imag->dx = stft->fRate;
    break;
  case RealStft:
  case PhaseStft:
  case HighResStft:
    stft->type = stftType;
    if(stft->real == NULL)
      stft->real = NewSignal();
    SizeSignal(stft->real,stft->nFrames*stft->nSubBands,YSIG); 
    ZeroSig(stft->real);
    stft->real->dx = stft->fRate;
    // We should also put the right firstp,lastp here
    // but this depends on depthHighRes !
    break;
  case HarmoStft:
    stft->type =  HarmoStft;
    
    // Compute the actual range of fundamental frequencies
    /* ***************************************************************/