stft.c

LastWave

    /*
     *                  Orthogonality condition : 
     *
     *  Given the freq0Id range [freq0IdMin,freq0IdMax], and the <windowSize>,
     *  what is the region of allowed indexes ON THE GRID i_f0=freq0Id/fRate,
     *  knowing that there is a LOWER BOUND i_f0 >= iFMO (iF0Min"Ortho")
     *  for these indexes in order to keep the partials at frequencies 
     *  k*freq0Id and (k+1)*freq0Id QUASI ORTHOGONAL.
     *  -> i_f0Min and i_f0Max depend on <windowSize> 
     */
    /* ***************************************************************/
    
    /* the f0 range */
    QuantizeRangeLarge(freq0IdMin,freq0IdMax,stft->fRate,
		       &f0IdMin,&f0IdMax);      
    i_f0Min = f0IdMin/stft->fRate;
    i_f0Max = f0IdMax/stft->fRate;
    
    i_f0Min = MAX(iFMO,i_f0Min);
    if(i_f0Min*stft->fRate < freq0IdMin) {
      i_f0Min++;
    }
    i_f0Max = MIN(stft->nSubBands-1,i_f0Max);
    
    /* This is now the freq0Id-range */
    stft->freqIdMin = i_f0Min*stft->fRate;
    stft->freqIdMax = i_f0Max*stft->fRate;
    stft->iFMO      = iFMO;
    stft->dimHarmo  = dimHarmo;
    
    /* If the range is not empty we allocate the data accordingly */
    if(stft->freqIdMin <= stft->freqIdMax) {
      stft->real = NewSignal();
      size = stft->nFrames*(((stft->freqIdMax-stft->freqIdMin)/stft->fRate)+1);
      SizeSignal(stft->real,size,YSIG); 
      ZeroSig(stft->real);
      stft->real->x0 = stft->freqIdMin;
      stft->real->dx = stft->fRate;
    }
    else {
      stft->real = NULL;
      Warningf("SizeStft : empty 'harmo' stft at windowSize=%d",windowSize);
    }
    
    break;
  default : 
    Errorf("SizeStft : unknown stftType %d",stftType);
  }
}   


/******************************************
 *
 *  Functions to access the real and imaginary
 *  inner product of a STFT structure
 *  corresponding to a given timeId (for all frequencies)
 *
 ******************************************/
void GetStftData(STFT stft,LWFLOAT timeId,LWFLOAT* *pReal,LWFLOAT* *pImag)
{
  int tId;
  
  /* Checking argument */
  CheckStft(stft);
  if(stft->real == NULL)
    Errorf("GetStftData : empty stft->real");
  if(pReal == NULL)
    Errorf("GetStftData : NULL output");
  if(pImag != NULL && stft->imag == NULL) 
    Errorf("GetStftData : trying to get imag in real stft");

  /* Some tests */
  if(timeId != (int) timeId) 
    Errorf("GetStftData : Bad timeId %g not integer",timeId);
  tId = (int) timeId;
  
  if (tId % stft->tRate != 0 || tId < 0 || tId/stft->tRate >= stft->nFrames) 
    Errorf("GetStftData : Bad timeId %d : rate %d range [0 %d]\n",
	   tId,stft->tRate,stft->nFrames);
  
  /* Getting the data */
  *pReal = stft->real->Y+(tId/stft->tRate)*(((stft->freqIdMax-stft->freqIdMin)/stft->fRate)+1);
  if(pImag != NULL)
    *pImag = stft->imag->Y+(tId/stft->tRate)*stft->nSubBands;
} 

/*
 * The GENERAL command for computing a stft
 */
 
void C_Stftd(char **argv)
{
  // The signal to be analyzed
  SIGNAL signal;
  struct aTFContent tfContent;
  // The output stft 
  STFT stft = NULL;

  // Some intermediate STFT variables
  STFT stftTmp        = NULL;
  STFT stftComplex    = NULL;
  STFT stftReal       = NULL;
  STFT stftPhase      = NULL;
  // Parameters for the grid

  int timeRedund,freqRedund;
  // The window
  char windowShape;
  int  windowSize;
  // Border effects
  char borderType;

  // The limits of the update
  int  timeIdMin,timeIdMax;
  
  char *name;
  char opt;
  char stftType;

  
  // Read the main arguments
  argv = ParseArgv(argv,tSTFT_,NULL,&stft,tSIGNALI,&signal,-1);
  if(stft == NULL) 
    stft = GetStftCur();
  // The signal size
  tfContent.signalSize = signal->size;
  tfContent.x0 = signal->x0;
  tfContent.dx = signal->dx;

  // The minimum and maximum timeId
  timeIdMin = 0;
  timeIdMax = tfContent.signalSize-1;
  
  // Default Option values
  windowShape = STFT_DEFAULT_WINDOWSHAPE;
  stftType    = STFT_DEFAULT_STFTTYPE;
  timeRedund  = STFT_DEFAULT_TIMEREDUND;
  freqRedund  = STFT_DEFAULT_FREQREDUND;
  borderType  = STFT_DEFAULT_BORDERTYPE;

  // Try to read the windowShape
  ParseArgv(argv,tSTR_,NULL,&name,-1);
  if(name != NULL) {
    windowShape = Name2WindowShape(name);
    argv++;
  }
  // The windowSize
  argv = ParseArgv(argv,tINT,&windowSize,-1);

  // Try to read the type
  ParseArgv(argv,tSTR_,NULL,&name,-1);
  if(name != NULL) {
    if(!strcmp(name,"complex") || !strcmp(name,"real")    || !strcmp(name,"phase")) {
      stftType = Name2StftType(name);
      argv++;
    }
  }

  // Reading additional options
  while(( opt = ParseOption(&argv))) {
    switch(opt) {
      // Border type
    case 'b' :
      argv = ParseArgv(argv,tSTR,&name,-1);
      borderType = Name2BorderType(name);
      break;
      // Time redundancy factor
    case 'T' :
      argv = ParseArgv(argv,tINT,&timeRedund,-1);
      if(timeRedund<=0) Errorf("Bad timeRedund %d",timeRedund);
      break;
      // Frequency redundancy factor
    case 'F' :
      argv = ParseArgv(argv,tINT,&freqRedund,-1);
      if(freqRedund<=0) Errorf("Bad freqRedund %d",freqRedund);
      break;
    default : ErrorOption(opt);
    }
  }
  NoMoreArgs(argv);

  //
  // Computation of a Complex spectrogram 
  //
  
  // -If we are computing a complex spectrogram 
  // we don't need any temporary stft structure
  if(stftType == ComplexStft) 
    stftTmp = stft;
  // -Otherwise we need one to store the (main) complex spectrogram.
  else
    stftTmp = NewStft();
  CopyFieldsTFContent(&tfContent,stftTmp);
  SizeStft(stftTmp,windowShape,windowSize,timeRedund,freqRedund,
	   borderType,ComplexStft,0,0,0,0);

  // Computing the complex spectrogram
  UpdateStftComplex(stftTmp,signal,timeIdMin,timeIdMax);
  // Are we done ?
  if(stftType == ComplexStft) {
    return;
  } else
    stftComplex = stftTmp;
  
  //
  // Computation of a Real spectrogram 
  //
  
  // -If we are computing a real or phase spectrogram
  // we don't need any temporary stft structure
  if(stftType == RealStft || stftType == PhaseStft) 
    stftTmp = stft;
  // -Otherwise we need one to store the real spectrogram.
  else
    stftTmp = NewStft();
  CopyFieldsTFContent(&tfContent,stftTmp);
  // If we compute ONLY the phase
  if(stftType == PhaseStft)
    SizeStft(stftTmp,windowShape,windowSize,timeRedund,freqRedund,borderType,
	     PhaseStft,0,0,0,0);
  // Else we need the RealStft in any case
  else
    SizeStft(stftTmp,windowShape,windowSize,timeRedund,freqRedund,borderType,
	     RealStft,0,0,0,0);

  // Compute the real or phase stft
  UpdateStftRealOrPhase(stftTmp,stftComplex,timeIdMin,timeIdMax); 
  // Are we done ?
  if (stftType == RealStft || stftType == PhaseStft) {
    // We delete the stftComplex 
    stftComplex = DeleteStft(stftComplex);
    return;
  }
  else
    stftReal = stftTmp;

}

/*
 * The TFContent fields : x0,dx,signalSize,freqIdNyquist
 */
static char *dxDoc = "{[= <dx>]} {Sets/Gets the abscissa step of the signal which has been analyzed.}";

void *GetDxTFContentV(ATFCONTENT val, void **arg)
{
  /* Documentation */
  if (val == NULL) return(dxDoc);
  
  return(GetFloatField(val->dx,arg));
}

void *SetDxTFContentV(ATFCONTENT val, void **arg)
{
  /* Documentation */
  if (val == NULL) return(dxDoc);

  return(SetFloatField(&(val->dx),arg,FieldSPositive));
}

static char *x0Doc = "{[= <x0>]} {Sets/Gets the first abscissa of the signal which has been analyzed.}";
void *GetX0TFContentV(ATFCONTENT val, void **arg)
{
  /* Documentation */
  if (val == NULL) return(x0Doc);
  return(GetFloatField(val->x0,arg));
}

void *SetX0TFContentV(ATFCONTENT val, void **arg)
{
  /* Documentation */
  if (val == NULL) return(x0Doc);
  return(SetFloatField(&(val->x0),arg,0));
}

static char *signalSizeDoc = "{} {Gets the size of the original signal of the time-frequency transform.}";
void *GetSignalSizeTFContentV(ATFCONTENT val, void **arg)
{
  /* Documentation */
  if (val == NULL) return(signalSizeDoc);
  return(GetIntField(val->signalSize,arg));
}

static char *freqIdNyquistDoc = "{} {Gets the Nyquist frequency in sample coordinates.}";

/* 
 * This is a bit special ... val is not used : 'freqIdNyquist' is a builtin constant 
 * that is useful for index<->real frequency conversions
*/
void *GetFreqIdNyquistTFContentV(VALUE val, void **arg)
{
  char *field = ARG_G_GetField(arg);
  /* Documentation */
  if (val == NULL)  return(freqIdNyquistDoc);
  return(GetIntField(GABOR_NYQUIST_FREQID,arg));
}

/*
 * The window fields : windowShape,windowSize
 */
static char *windowShapeDoc = "{} {Gets the windowShape of a Short Time Fourier Transform. "
WindowShapeHelpString
". To see the shape of a window, you can build it using the 'stft window ...' command.}";
static char *windowSizeDoc = "{} {Gets the windowSize of a Short Time Fourier Transform, i.e. the number of samples of its window. So far only powers of 2 are allowed.  To see the shape of a window, you can build it using the 'stft window ...' command.}";

void *GetWindowShapeStftV(STFT val, void **arg)
{
  /* Documentation */
  if (val == NULL) return(windowShapeDoc);
  return(GetStrField(WindowShape2Name(val->windowShape),arg));
}

void *GetWindowSizeStftV(STFT val, void **arg)
{
  /* Documentation */
  if (val == NULL) return(windowSizeDoc);
  return(GetIntField(val->windowSize,arg));
}

static char *gridDoc = "{} {Gets stft 'grid' i.e. returns a &listv {timeRate timeLength freqRate freqLength}.}";
void *GetGridStftV(STFT stft,void **arg)
{
  LISTV lv;
  if(stft == NULL) return(gridDoc);

  lv = TNewListv();
  AppendInt2Listv(lv,stft->tRate);
  AppendInt2Listv(lv,stft->nFrames);
  AppendInt2Listv(lv,stft->fRate);
  AppendInt2Listv(lv,stft->nSubBands);
  return(GetValueField(lv,arg));
}

static char *typeDoc = "{} {Gets stft type (it is either 'complex' for short time fourier transform, 'real', 'phase' or 'highres'}";
void *GetTypeStftV(STFT val,void **arg)
{
  if(val == NULL) return(typeDoc);
  return(GetStrField(StftType2Name(val->type),arg));
}

static char *borderDoc = "{} {Returns stft <border type>"
BorderTypeHelpString
".}";
static char *firstpDoc = "{} {Gets stft <firstp>}";
static char *lastpDoc = "{} {Gets stft <lastp>}";

void *GetBorderStftV(STFT stft,void **arg)
{
  if(stft == NULL) return(borderDoc);
  return(GetStrField(BorderType2Name(stft->borderType),arg));
}

void *GetFirstpStftV(STFT stft,void **arg)
{
  if(stft == NULL) return(firstpDoc);
  return(GetIntField(stft->firstp,arg));
}

void *GetLastpStftV(STFT stft,void **arg)
{
  if(stft == NULL) return(lastpDoc);
  return(GetIntField(stft->lastp,arg));
}
  
static char *dataDoc = "{} {Returns a signal containing the stft real data or a listv {real imag} containing its complex data}";
void *GetDataStftV(STFT stft,void **arg)
{
  SIGNAL signalR = NULL,signalI = NULL;
  LISTV lv = NULL;

  /* Documentation */
  if(stft == NULL) return(dataDoc);
  
  /* Case of an empty stft */
  if(stft->real==NULL) {
    return(GetValueField(nullValue,arg));
  }
  switch(stft->type) {
  case ComplexStft :
    lv = TNewListv();
    signalR = TNewSignal();
    signalI = TNewSignal();
    CopySig(stft->real,signalR);
    CopySig(stft->imag,signalI);
    signalR->x0 = FreqId2Freq(stft,stft->freqIdMin);
    signalR->dx = FreqId2Freq(stft,stft->fRate);
    signalI->x0 = FreqId2Freq(stft,stft->freqIdMin);
    signalI->dx = FreqId2Freq(stft,stft->fRate);
    AppendValue2Listv(lv,(VALUE)signalR);
    AppendValue2Listv(lv,(VALUE)signalI);
    return(GetValueField(lv,arg));
  case RealStft: 
  case PhaseStft: 
  case HighResStft:
  case HarmoStft:
    signalR = TNewSignal();
    CopySig(stft->real,signalR);
    signalR->x0 = FreqId2Freq(stft,stft->freqIdMin);
    signalR->dx = FreqId2Freq(stft,stft->fRate);
    return(GetValueField(signalR,arg));
  }
}

/*
 * The field list
 */
struct field fieldsStft[] = {
  "dx",GetDxTFContentV,SetDxTFContentV,NULL,NULL,
  "x0",GetX0TFContentV,SetX0TFContentV,NULL,NULL,
  "signalSize",GetSignalSizeTFContentV,NULL,NULL,NULL,
  "freqIdNyquist",GetFreqIdNyquistTFContentV,NULL,NULL,NULL,
  "windowShape",GetWindowShapeStftV,NULL,NULL,NULL,
  "windowSize",GetWindowSizeStftV,NULL,NULL,NULL,
  "grid",GetGridStftV,NULL,NULL,NULL,
  "border",GetBorderStftV,NULL,NULL,NULL,
  "firstp",GetFirstpStftV,NULL,NULL,NULL,
  "lastp",GetLastpStftV,NULL,NULL,NULL,
  "type",GetTypeStftV,NULL,NULL,NULL,
  "sig",GetDataStftV,NULL,NULL,NULL,
  NULL, NULL, NULL, NULL, NULL
};

/*
 * The type structure for STFT
 */

TypeStruct tsStft = {

  "{{{&stft} {This type is the basic type for Short Time Fourier transforms and related time-frequency transforms.}}}",  /* Documentation */

  &stftType,       /* The basic (unique) type name */
  NULL,     /* The GetType function */                       
  
  DeleteStft,     /* The Delete function */
  NewStft,     /* The New function */
  
  NULL,       /* The copy function */
  ClearStft,       /* The clear function */
  
  ToStrStft,       /* String conversion */
  ShortPrintStft,   /* The Print function : print the object when 'print' is called */
  PrintInfoStft,   /* The PrintInfo function : called by 'info' */

  NumExtractStft,              /* The NumExtract function : used to deal with syntax like 10a */
   
  fieldsStft,      /* The list of fields */
};

/* EOF */