stft.c

LastWave

void StftOper(const STFT stftIn1,const STFT stftIn2,STFT stftOut,char oper)
{
  int i;
  LWFLOAT re1,im1,re2,im2,norm2;
  static double log2 = -1;

  if(log2 == -1) log2 = log(2);

  CheckStftNotEmpty(stftIn1);
  if(stftIn2 != NULL) {
    CheckStftNotEmpty(stftIn2);
    CheckSameStftStruct(stftIn1,stftIn2);
  }
  CheckStft(stftOut);

  if(!StftOperIsOK(oper))
    Errorf("StftOper : unknown arithmetic operation");

  // If the output is different from the inputs we size it
  if(stftOut != stftIn1 && stftOut != stftIn2) {
    CopyFieldsTFContent(stftIn1,stftOut);
    // TODO : deal with time/freqRedund ... should be computed from tRate ?
    SizeStft(stftOut,stftIn1->windowShape,stftIn1->windowSize,
	     STFT_DEFAULT_TIMEREDUND,STFT_DEFAULT_FREQREDUND,
	     stftIn1->borderType,stftIn1->type,
	     stftIn1->freqIdMin,stftIn1->freqIdMax,
	     stftIn1->iFMO,stftIn1->dimHarmo);
  }
  // DEBUG
  CheckSameStftStruct(stftIn1,stftOut);

  switch(stftIn1->type) {
  case ComplexStft:
    switch(oper) {
    case STFT_ADD: //Should work even if stftOut==stftIn1 or stftIn2
      CheckStftNotEmpty(stftIn2);
      for(i = 0; i < stftIn1->real->size; i++) {
	stftOut->real->Y[i] = stftIn1->real->Y[i]+stftIn2->real->Y[i];
	stftOut->imag->Y[i] = stftIn1->imag->Y[i]+stftIn2->imag->Y[i];
      }
      break;
    case STFT_SUBSTRACT: //Should work even if stftOut==stftIn1 or stftIn2
      CheckStftNotEmpty(stftIn2);
      for(i = 0; i < stftIn1->real->size; i++) {
	stftOut->real->Y[i] = stftIn1->real->Y[i]-stftIn2->real->Y[i];
	stftOut->imag->Y[i] = stftIn1->imag->Y[i]-stftIn2->imag->Y[i];
      }
      break;
    case STFT_MULTIPLY: //Should work even if stftOut==stftIn1 or stftIn2
      CheckStftNotEmpty(stftIn2);
      for(i = 0; i < stftIn1->real->size; i++) {
	re1 = stftIn1->real->Y[i];
	re2 = stftIn2->real->Y[i];
	im1 = stftIn1->imag->Y[i];
	im2 = stftIn2->imag->Y[i];
	stftOut->real->Y[i] = re1*re2-im1*im2;
	stftOut->imag->Y[i] = re1*im2+im1*re2;
      }
      break;
    case STFT_DIVIDE: //Should work even if stftOut==stftIn1 or stftIn2
      CheckStftNotEmpty(stftIn2);
      for(i = 0; i < stftIn1->real->size; i++) {
	re1 = stftIn1->real->Y[i];
	re2 = stftIn2->real->Y[i];
	im1 = stftIn1->imag->Y[i];
	im2 = stftIn2->imag->Y[i];
	norm2 = re2*re2+im2*im2;
	if(norm2 == 0) {
	  Warningf("StftOper : division by zero");
	  stftOut->real->Y[i] = FLT_MAX/2;
	  stftOut->imag->Y[i] = FLT_MAX/2;
	}
	else {
	  stftOut->real->Y[i] = (re1*re2+im1*im2)/norm2;
	  stftOut->imag->Y[i] = (-re1*im2+im1*re2)/norm2;
	}
      }
      break;
    case STFT_CONJUGATE:
      if(stftIn2 != NULL)
	Errorf("StftOper : stftIn2 should be NULL for '%s' operation",StftOper2Name(oper));
      for(i = 0; i < stftIn1->real->size; i++) {
	stftOut->real->Y[i] = stftIn1->real->Y[i];
	stftOut->imag->Y[i] = -stftIn1->imag->Y[i];
      }
      break;
    case STFT_LN: //Natural logarithm
      if(stftIn2 != NULL)
	Errorf("StftOper : stftIn2 should be NULL for '%s' operation",StftOper2Name(oper));
      for(i = 0; i < stftIn1->real->size; i++) {
	re1 = stftIn1->real->Y[i];
	im1 = stftIn1->imag->Y[i];
	norm2 = re1*re1+im1*im1;
	if(norm2 == 0.0) {
	  Warningf("StftOper : log of zero");
	  stftOut->real->Y[i]  = -FLT_MAX/2;
	  stftOut->imag->Y[i] = 0.0;
	}
	else {
	  stftOut->real->Y[i] = 0.5*log(norm2);
	  stftOut->imag->Y[i] = atan2(im1,re1)/(2*M_PI);
	}
      }
      break;
    case STFT_LOG: //Base 10
      if(stftIn2 != NULL)
	Errorf("StftOper : stftIn2 should be NULL for '%s' operation",StftOper2Name(oper));
      for(i = 0; i < stftIn1->real->size; i++) {
	re1 = stftIn1->real->Y[i];
	im1 = stftIn1->imag->Y[i];
	norm2 = re1*re1+im1*im1;
	if(norm2 == 0.0) {
	  Warningf("StftOper : log10 of zero");
	  stftOut->real->Y[i]  = -FLT_MAX/2;
	  stftOut->imag->Y[i] = 0.0;
	}
	else {
	  stftOut->real->Y[i] = 0.5*log10(norm2);
	  stftOut->imag->Y[i] = atan2(im1,re1)/(2*M_PI);
	}
      }
      break;
    case STFT_LOG2: //Base 2
      if(stftIn2 != NULL)
	Errorf("StftOper : stftIn2 should be NULL for '%s' operation",StftOper2Name(oper));
      for(i = 0; i < stftIn1->real->size; i++) {
	re1 = stftIn1->real->Y[i];
	im1 = stftIn1->imag->Y[i];
	norm2 = re1*re1+im1*im1;
	if(norm2 == 0.0) {
	  Warningf("StftOper : log2 of zero");
	  stftOut->real->Y[i]  = -FLT_MAX/2;
	  stftOut->imag->Y[i] = 0.0;
	}
	else {
	  stftOut->real->Y[i] = 0.5*log(norm2)/log2;
	  stftOut->imag->Y[i] = atan2(im1,re1)/(2*M_PI);
	}
      }
      break;
    default :
      Errorf("StftOper : operation '%s' not treated fro '%s' stft",StftOper2Name(oper),StftType2Name(stftIn1->type));
      break;
    }
    break;
  case RealStft:
  case PhaseStft:
  case HighResStft:
  case HarmoStft:
    switch(oper) {
    case STFT_ADD://Should work even if stftOut==stftIn1 or stftIn2
      CheckStftNotEmpty(stftIn2);
      for(i = 0; i < stftIn1->real->size; i++) {
	stftOut->real->Y[i] = stftIn1->real->Y[i]+stftIn2->real->Y[i];
      }
      break;
    case STFT_SUBSTRACT://Should work even if stftOut==stftIn1 or stftIn2
      CheckStftNotEmpty(stftIn2);
      for(i = 0; i < stftIn1->real->size; i++) {
	stftOut->real->Y[i] = stftIn1->real->Y[i]-stftIn2->real->Y[i];
      }
      break;
    case STFT_MULTIPLY://Should work even if stftOut==stftIn1 or stftIn2
      CheckStftNotEmpty(stftIn2);
      for(i = 0; i < stftIn1->real->size; i++) {
	stftOut->real->Y[i] = stftIn1->real->Y[i]*stftIn2->real->Y[i];
      }
      break;
    case STFT_DIVIDE://Should work even if stftOut==stftIn1 or stftIn2
      CheckStftNotEmpty(stftIn2);
      for(i = 0; i < stftIn1->real->size; i++) {
	if(stftIn2->real->Y[i] == 0) {
	  Warningf("StftOper : division by zero");
	  stftOut->real->Y[i] = FLT_MAX/2;
	}
	else {
	  stftOut->real->Y[i] = stftIn1->real->Y[i]/stftIn2->real->Y[i];
	}
      }
      break;
    case STFT_LN: //Natural logarithm
      if(stftIn2 != NULL)
	Errorf("StftOper : stftIn2 should be NULL for '%s' operation",StftOper2Name(oper));
      for(i = 0; i < stftIn1->real->size; i++) {
	re1 = stftIn1->real->Y[i];
	if(re1 <= 0.0) {
	  Warningf("StftOper : log of %g",re1);
	  stftOut->real->Y[i]  = -FLT_MAX/2;
	}
	else {
	  stftOut->real->Y[i] = log(re1);
	}
      }
      break;
    case STFT_LOG: //Base 10 logarithm
      if(stftIn2 != NULL)
	Errorf("StftOper : stftIn2 should be NULL for '%s' operation",StftOper2Name(oper));
      for(i = 0; i < stftIn1->real->size; i++) {
	re1 = stftIn1->real->Y[i];
	if(re1 <= 0.0) {
	  Warningf("StftOper : log10 of %g",re1);
	  stftOut->real->Y[i]  = -FLT_MAX/2;
	}
	else {
	  stftOut->real->Y[i] = log10(re1);
	}
      }
      break;
    case STFT_LOG2: //Base 2 logarithm
      if(stftIn2 != NULL)
	Errorf("StftOper : stftIn2 should be NULL for '%s' operation",StftOper2Name(oper));
      for(i = 0; i < stftIn1->real->size; i++) {
	re1 = stftIn1->real->Y[i];
	if(re1 <= 0.0) {
	  Warningf("StftOper : log2 of %g",re1);
	  stftOut->real->Y[i]  = -FLT_MAX/2;
	}
	else {
	  stftOut->real->Y[i] = log(re1)/log(2);
	}
      }
      break;
    default :
      Errorf("StftOper : operation '%s' not treated for '%s' stft",StftOper2Name(oper),StftType2Name(stftIn1->type));
      break;
    }
    break;
  default : 
    Errorf("StftOper : stft type '%s' not treated yet",StftType2Name(stftIn1->type));
  }
  stftOut->flagUpToDate = YES;
}


/* COMMAND */
void C_Stft(char **argv)
{
  STFT stft = NULL;
  STFT stft1 = NULL;
  STFT stft2 = NULL;
  STREAM stream;
  char *action;
  char opt;
  
  char flagShortPrint = NO;
  char flagHeader = YES;
  char *filename;
  
  // To get tabulated data from the package
  SIGNAL signal,signal1,real,imag;
  char *windowShapeName;
  unsigned long windowSize;
  LISTV listv = NULL;

  /* Parsing command */
  argv = ParseArgv(argv,tWORD,&action,-1);
  
  if(!strcmp(action,"window")) {
    argv = ParseArgv(argv,tSTR,&windowShapeName,tINT,&windowSize,0);
    GetTabWindow(Name2WindowShape(windowShapeName),windowSize,&real);
    signal = TNewSignal();
    CopySig(real,signal);
    SetResultValue(signal);
    return;
  }
  if(!strcmp(action,"exp")) {
    NoMoreArgs(argv);
    GetTabExponentials(&real,&imag);
    signal= TNewSignal();
    signal1= TNewSignal();
    CopySig(real,signal);
    CopySig(imag,signal1);
    listv = TNewListv();
    AppendValue2Listv(listv,(VALUE)signal);
    AppendValue2Listv(listv,(VALUE)signal1);
    SetResultValue(listv);
    return;
  }
  if(!strcmp(action,"gg")) {
    argv = ParseArgv(argv,tSTR,&windowShapeName,tINT,&windowSize,0);
    GetTabGG(Name2WindowShape(windowShapeName),windowSize,&real,&imag);
    signal= TNewSignal();
    signal1= TNewSignal();
    CopySig(real,signal);
    CopySig(imag,signal1);
    listv = TNewListv();
    AppendValue2Listv(listv,(VALUE)signal);
    AppendValue2Listv(listv,(VALUE)signal1);
    SetResultValue(listv);
    return;
  }

  if(!strcmp(action,"write")) {
    argv = ParseArgv(argv,tSTFT,&stft,-1);
    argv = ParseArgv(argv,tSTREAM_,NULL,&stream,-1);
    if (stream == NULL) argv = ParseArgv(argv,tSTR,&filename,-1);
    else filename = NULL;
    while((opt = ParseOption(&argv))) { 
      switch(opt) {
      case 'h': flagHeader = NO; break;
      default: ErrorOption(opt);
      }
    }    
    NoMoreArgs(argv);
    if (stream == NULL) WriteStftFile(stft,filename,flagHeader);
    else WriteStftStream(stft,stream,flagHeader);
    return;
  }
  
  if(!strcmp(action,"+") || !strcmp(action,"-") || 
     !strcmp(action,"*") || !strcmp(action,"/")) {
    argv = ParseArgv(argv,tSTFT,&stft,tSTFT,&stft1,tSTFT_,NULL,&stft2,0);
    // If there are only two arguments
    if (stft2 == NULL) stft2 = stft;
    StftOper(stft,stft1,stft2,Name2StftOper(action));
    return;
  }

  if(!strcmp(action,"conjugate") || !strcmp(action,"ln") ||
     !strcmp(action,"log") || !strcmp(action,"log2")) {
    argv = ParseArgv(argv,tSTFT,&stft,tSTFT_,NULL,&stft1,0);
    // If there is only one argument
    if (stft1 == NULL) stft1 = stft;
    StftOper(stft,NULL,stft1,Name2StftOper(action));
    return;
  }
  Printf("Unknown action '%s'\n",action);
  ErrorUsage();
}

/************************************/
/*
 *	Type conversions 
 */				    
/************************************/
char * StftType2Name(char stftType)
{
  if(!StftTypeIsOK(stftType))
    Errorf("StftType2Name : bad stftType %d",stftType);
  
  switch(stftType) {
  case ComplexStft: return("complex");
  case RealStft:    return("real");
  case PhaseStft:   return("phase");
  case HighResStft: return("highres");
  case HarmoStft:   return("harmo");
  default : Errorf("StftType2Name : (Weird) bad type %d",stftType);
  }
  return("");
}


char Name2StftType(char *stftType)
{
  if (!strcmp(stftType,"complex")) return(ComplexStft);
  if (!strcmp(stftType,"real")) return(RealStft);
  if (!strcmp(stftType,"phase")) return(PhaseStft);
  if (!strcmp(stftType,"highres")) return(HighResStft);
  if (!strcmp(stftType,"harmo")) return(HarmoStft);
  Errorf("Name2StftType : (Weird) bad stft name %s",stftType);
  // This should never be reached but the compiler may want to be sure the fnuction returns a value
  return(LastStftType);
}


/************************************/
/*
 *	ALLOCATIONS
 *
 */				    
/************************************/
static void _MyCheckHarmo(const STFT stft)
{
  if(!INRANGE(STFT_MIN_WINDOWSIZE,stft->windowSize,stft->signalSize))
    Errorf("CheckHarmo : bad windowSize %d not in [4 %d]",stft->windowSize,stft->signalSize);
  
  /* Checking frequency grid */
  if(stft->fRate*(stft->nSubBands-1) != GABOR_NYQUIST_FREQID)
    Errorf("CheckHarmo : bad freq grid %d * %d != %d for windowSize %d",
	   stft->fRate,stft->nSubBands-1,GABOR_NYQUIST_FREQID,stft->windowSize);
  
  /* The freq0Id-range must be on the grid */
  if(stft->freqIdMin % stft->fRate != 0 || stft->freqIdMax % stft->fRate != 0)
    Errorf("CheckHarmo : bad f0 range [%d %d] not on the grid %d",stft->freqIdMin,stft->freqIdMax,stft->fRate);
  
  if(stft->freqIdMin < stft->iFMO*stft->fRate)
    Errorf("CheckHarmo ; bad freq0IdMin %d < %d",stft->freqIdMin,stft->iFMO*stft->fRate);
  if(stft->freqIdMin <= stft->freqIdMax && stft->real == NULL)
    Errorf("CheckHarmo : NULL data for range [%d %d]",stft->freqIdMin,stft->freqIdMax);
  if(stft->freqIdMin > stft->freqIdMax && stft->real != NULL)
    Errorf("CheckHarmo : NON NULL data for range [%d %d]",stft->freqIdMin,stft->freqIdMax);
  if(stft->freqIdMax > stft->fRate*(stft->nSubBands-1))
    Errorf("CheckHarmo ; bad freq0IdMax %d > %d",stft->freqIdMax,stft->fRate*(stft->nSubBands-1));
}

void CheckStft(const STFT stft)
{
  // Checking of fields independently of one another
  if (!BorderTypeIsOK(stft->borderType))
    Errorf("CheckStft : unknown border type %d",stft->borderType); 
  if (!WindowShapeIsOK(stft->windowShape)) 
    Errorf("CheckStft : unknown windowShape %d",stft->windowShape); 
  if(stft->tRate <= 0 || stft->fRate <= 0) 
    Errorf("CheckStft : bad rates %d %d",stft->tRate,stft->fRate); 
  if(stft->tRate*stft->nFrames < stft->signalSize) 
    Errorf("CheckStft : bad time grid %d * %d < %d",stft->tRate,stft->nFrames,stft->signalSize);

  if(stft->iFMO <= 0)
    Errorf("CheckStft : bad iFMO %d for HarmoStft",stft->iFMO);
  if(stft->dimHarmo <= 0)
    Errorf("CheckStft : bad dimHarmo %d for HarmoStft",stft->dimHarmo);
  
  // Checking of coherence of fields
  /* Case of an empty stft */
  if(stft->signalSize == 0) {
    if(stft->windowSize > 0)
      Errorf("CheckStft %d %s : empty stft is at windowSize %d!",
	     stft->windowSize,StftType2Name(stft->type),stft->windowSize);
    if(stft->real != NULL || stft->imag != NULL)
      Errorf("CheckStft %d %s : non NULL data for empty stft",
	     stft->windowSize,StftType2Name(stft->type));
    if(stft->flagUpToDate)
      Errorf("CheckStft %d %s : empty stft is up to date!");
  }
  /* Case of non empty one */
  else {
    // Temporary during development
    if(stft->type == HarmoStft) {    
      _MyCheckHarmo(stft);
      return;
    }
    if(stft->windowSize <= 0)
      Errorf("CheckStft %d %s : bad windowSize %d <= 0",
	     stft->windowSize,StftType2Name(stft->type),stft->windowSize);
    /* Window size must be smaller than signal size ! */
    if(!INRANGE(2,stft->windowSize,stft->signalSize))
      Errorf("CheckStft %d %s : bad windowSize %d compared to signalSize %d",
	     stft->windowSize,StftType2Name(stft->type),stft->windowSize,stft->signalSize);
    /* TODO ; remove ??? */
    if(stft->windowSize > GABOR_MAX_FREQID)
      Errorf("CheckStft %d %s : bad windowSize %d compared to GABOR_MAX_FREQID %d",
	     stft->windowSize,StftType2Name(stft->type),stft->windowSize,GABOR_MAX_FREQID);
  }
  
  /* Checking types */
  if(!StftTypeIsOK(stft->type))
    Errorf("CheckStft %d %s : unknown  type %d",stft->windowSize,StftType2Name(stft->type),stft->type);