stft.c

LastWave

/*..........................................................................*/
/*                                                                          */
/*      L a s t W a v e    P a c k a g e 'stft' 2.1                         */
/*                                                                          */
/*      Copyright (C) 1997-2002 R.Gribonval, E.Bacry                        */
/*      email  : [e-mail]                                    */
/*               [e-mail]                            */
/*                                                                          */
/*..........................................................................*/
/*                                                                          */
/*      This program is a free software, you can redistribute it and/or     */
/*      modify it under the terms of the GNU General Public License as      */
/*      published by the Free Software Foundation; either version 2 of the  */
/*      License, or (at your option) any later version                      */
/*                                                                          */
/*      This program is distributed in the hope that it will be useful,     */
/*      but WITHOUT ANY WARRANTY; without even the implied warranty of      */
/*      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       */
/*      GNU General Public License for more details.                        */
/*                                                                          */
/*      You should have received a copy of the GNU General Public License   */
/*      along with this program (in a file named COPYRIGHT);                */
/*      if not, write to the Free Software Foundation, Inc.,                */
/*      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA             */
/*                                                                          */
/*..........................................................................*/



/****************************************************/
/*
 * 	The Short Time Fourier Transform :
 *		Parsing of variables
 *		Allocation and Desallocations of structures
 *
 */
/****************************************************/
#include "lastwave.h"

#include "stft.h"

/********************************/
/*
 *	STFT variables
 */
/********************************/
char *stftType = "&stft";
extern TypeStruct tsStft;


/*
 * Answers to the different print messages
 */
 
void ShortPrintStft(STFT stft)
{
  //  if (stft->name == NULL)
    Printf("<&stft[%d,%d];%p>\n",stft->nFrames,stft->nSubBands,stft);
    //  else
    //    Printf("<'%s';&stft[%d,%d];%p>\n",stft->name,stft->nFrames,stft->nSubBands,stft);
}

char *ToStrStft(STFT stft, char flagShort)
{
  static char str[30];
  
  sprintf(str,"<&stft;%p>",stft);
  return(str);
}


void PrintInfoStft(STFT stft)
{
  PrintStft(stft,NO);
}

/*
 * NumExtraction
 *
 * Vertical slices (<timeId>stft, ...)
 * Horizontal slices (<freqId>stft ...)
 */
static char *numdoc = "The syntax <timeId>stft corresponds to vertical slices of the &stft. The syntax .<freqId>stft addresses horizontal slices. A slice is either a &listv {<real> <imag>} if the type of the &stft is complex, or just a &signal <real> which may contain either the energy, the phase, or the highres correlation.";

static void *NumExtractStft(STFT stft,void **arg)
{
  long n;
  char flagDot;
  
  int timeId,freqId;
  SIGNAL signalR,signalI;
  LWFLOAT *real,*imag;
  LISTV lv = NULL;

  /* doc */
  if (stft == NULL) return(numdoc);

  n = ARG_NE_GetN(arg);
  flagDot = ARG_NE_GetFlagDot(arg);

  /* Horizontal slice at as fixed frequency */
  if(flagDot) {
    CheckStftNotEmpty(stft);
    freqId = n;
    /* Checking */
    if(freqId%stft->fRate != 0 || 
       !INRANGE(0,freqId/stft->fRate,stft->nSubBands-1) ||
       !INRANGE(stft->freqIdMin,freqId,stft->freqIdMax)) {
	SetErrorf("Bad freqId %d : rate %d length %d min %d max %d",freqId,stft->fRate,stft->nSubBands,stft->freqIdMin,stft->freqIdMax);
	return(NULL);
    }
      /* Checking */
    switch(stft->type) {
    case ComplexStft:
      signalR = TNewSignal();
      signalI = TNewSignal();
      SizeSignal(signalR,stft->nFrames,YSIG);
      SizeSignal(signalI,stft->nFrames,YSIG);
      signalR->x0      =signalI->x0      =stft->x0;
      signalR->dx      =signalI->dx      = TimeId2Time(stft,stft->tRate);
      signalR->firstp  =signalI->firstp  = stft->firstp;
      signalR->lastp   =signalI->lastp   = stft->lastp;
      /* Getting the data */
      for(timeId = 0;
	  timeId < stft->tRate*stft->nFrames;
	  timeId += stft->tRate) {
	GetStftData(stft,timeId,&real,&imag);
	signalR->Y[timeId/stft->tRate] = real[(freqId-stft->freqIdMin)/stft->fRate];
	signalI->Y[timeId/stft->tRate] = imag[(freqId-stft->freqIdMin)/stft->fRate];
      }
      lv = TNewListv();
      AppendValue2Listv(lv,(VALUE)signalR);
      AppendValue2Listv(lv,(VALUE)signalI);
      ARG_NE_SetResValue(arg,lv);
      return(listvType);
    case RealStft: 
    case PhaseStft: 
    case HighResStft:
    case HarmoStft:
      signalR = TNewSignal();
      SizeSignal(signalR,stft->nFrames,YSIG);
      signalR->x0 = stft->x0;
      signalR->dx = TimeId2Time(stft,stft->tRate);
      signalR->firstp  = stft->firstp;
      signalR->lastp   = stft->lastp;
      /* Getting the data */
      for(timeId = 0; 
	  timeId < stft->tRate*stft->nFrames; 
	  timeId += stft->tRate) {
	GetStftData(stft,timeId,&real,NULL);
	signalR->Y[timeId/stft->tRate] = real[(freqId-stft->freqIdMin)/stft->fRate];
      }
      ARG_NE_SetResValue(arg,signalR);
      return(signalType);
    }
  } else {
    /* Vertical slice at as fixed time */
    CheckStftNotEmpty(stft);
    timeId = n;
    switch(stft->type) {
    case ComplexStft:
      /* Getting the data */
      GetStftData(stft,timeId,&real,&imag);
      /* Allocation */
      signalR = TNewSignal();
      signalI = TNewSignal();
      SizeSignal(signalR,stft->nSubBands,YSIG);
      SizeSignal(signalI,stft->nSubBands,YSIG);
      signalR->x0 = signalI->x0 = 0.0;
      signalR->dx = signalI->dx = FreqId2Freq(stft,stft->fRate);
      signalR->firstp  = signalI->firstp  = 0;
      signalR->lastp   = signalI->lastp   = stft->nSubBands-1;
      /* Copying */
      memcpy(signalR->Y,real,stft->nSubBands*sizeof(LWFLOAT));
      memcpy(signalI->Y,imag,stft->nSubBands*sizeof(LWFLOAT));
      lv = TNewListv();
      AppendValue2Listv(lv,(VALUE)signalR);
      AppendValue2Listv(lv,(VALUE)signalI);
      ARG_NE_SetResValue(arg,lv);
      return(listvType);
    case RealStft:
    case PhaseStft:
    case HighResStft:
    case HarmoStft:
      /* Getting the data */
      GetStftData(stft,timeId,&real,NULL);
      /* Allocation */
      signalR = TNewSignal();
      SizeSignal(signalR,((stft->freqIdMax-stft->freqIdMin)/stft->fRate)+1,YSIG);
      signalR->x0      = FreqId2Freq(stft,stft->freqIdMin);
      signalR->dx      = FreqId2Freq(stft,stft->fRate);
      signalR->firstp  = 0;
      signalR->lastp   = stft->nSubBands-1;
      /* Copying */
      memcpy(signalR->Y,real,(((stft->freqIdMax-stft->freqIdMin)/stft->fRate)+1)*sizeof(LWFLOAT));
      ARG_NE_SetResValue(arg,signalR);
      return(signalType);
    }
  }
}



STFT TNewStft(void)
{
  STFT stft;
  
  stft = NewStft();
  TempValue(stft);
  return(stft);
}

/*
 * Get the current stft
 * (generate an error if there is none)
 */

STFT GetStftCur(void)
{
  STFT stft;

  if(!ParseTypedValLevel_(levelCur,"objCur",NULL,(VALUE *)&stft,stftType)) 
    Errorf1("");
  
  AddRefValue(stft);
  TempValue(stft);
  
  return(stft);
}


static void CheckSameStftStruct(const STFT stft1,const STFT stft2)
{
  if(stft1==NULL||stft2==NULL)
    Errorf("CheckSameStftStruct : NULL input");
  CheckTFContentCompat(stft1,stft2);
  if(stft1->tRate != stft2->tRate ||
     stft1->nFrames != stft2->nFrames ||
     stft1->fRate != stft2->fRate ||
     stft1->nSubBands != stft2->nSubBands ||
     stft1->borderType != stft2->borderType ||
     stft1->firstp != stft2->firstp ||
     stft1->lastp != stft2->lastp ||
     stft1->freqIdMin != stft2->freqIdMin ||
     stft1->freqIdMax != stft2->freqIdMax ||
     stft1->windowShape != stft2->windowShape ||
     stft1->windowSize != stft2->windowSize ||
     stft1->dimHarmo != stft2->dimHarmo ||
     stft1->iFMO != stft2->iFMO ||
     stft1->type != stft2->type) {
    PrintStft(stft1,NO);
    PrintStft(stft2,NO);
    Errorf("CheckSameStftStruct : struct is different!");
  }
}

char* StftOper2Name(char oper)
{
  switch(oper) {
  case STFT_ADD:
    return("+");
  case STFT_SUBSTRACT:
    return("-");
  case STFT_MULTIPLY:
    return("*");
  case STFT_DIVIDE:
    return("/");
  case STFT_CONJUGATE:
    return("conjugate");
  case STFT_LN:
    return("ln");
  case STFT_LOG:
    return("log");
  case STFT_LOG2:
    return("log2");
  default :
    Errorf("StftOper2Name : (Weird) unknown oper %s",oper);
    return("");
  }
}

char Name2StftOper(char *name)
{
  if(name==NULL)
    Errorf("Name2StftOper : NULL input");
  if(!strcmp(name,"+"))
    return(STFT_ADD);
  if(!strcmp(name,"-"))
    return(STFT_SUBSTRACT);
  if(!strcmp(name,"*"))
    return(STFT_MULTIPLY);
  if(!strcmp(name,"/"))
    return(STFT_DIVIDE);
  if(!strcmp(name,"conjugate"))
    return(STFT_CONJUGATE);
  if(!strcmp(name,"ln"))
    return(STFT_LN);
  if(!strcmp(name,"log"))
    return(STFT_LOG);
  if(!strcmp(name,"log2"))
    return(STFT_LOG2);
  Errorf("Name2StftOper : unknown operation '%s'",name);
  return(-1);
}


// Given a stft and a time range [timeIdMin,timeIdMax] where a
// signal has changed, computes the time range where to update the stft
void  ComputeUpdateLimits(STFT stft,int timeIdMin,int timeIdMax,
			  int *pTimeIdMin,int *pTimeIdMax,
			  char *pFlagTwoLoops,int *pTimeIdMin1,int *pTimeIdMax1)
{
  int shiftMin,shiftMax;
  // Check ...
  CheckStftNotEmpty(stft);
  switch(stft->borderType) {
    // The simple border effects
  case BorderPad0 :
  case BorderPad : 
  case BorderMirror :
    if((!INRANGE(0,timeIdMin,timeIdMax)) || 
       (!INRANGE(timeIdMin,timeIdMax,stft->signalSize-1)))
      Errorf("ComputeUpdateLimits : [%d %d] not in [0 %d]",
	     timeIdMin,timeIdMax,stft->signalSize-1);
    break;
  case BorderPeriodic :
    if((!INRANGE(0,timeIdMin,MIN(timeIdMax,stft->signalSize-1))) || 
       (!INRANGE(timeIdMin,timeIdMax,2*stft->signalSize-1)))
      Errorf("ComputeUpdateLimits : [%d %d] not in [0 %d]",
	     timeIdMin,timeIdMax,stft->signalSize-1);
    break;
  default : 
    Errorf("ComputeUpdateLimits : border '%s' not treated yet",BorderType2Name(stft->borderType));
  }

  // By default there is only one loop
  *pFlagTwoLoops = NO;

  // The time support of the window centered at timeId is
  // [timeId+shiftMin,timeId+shiftMax]
  ComputeWindowSupport(stft->windowSize,stft->windowShape,0.0,
		       &shiftMin,&shiftMax);
  // The inner product has changed iff this support intersects 
  // the portion of the signal that has changed. This depends on the border type
  switch(stft->borderType) {
    // The simple border effects
  case BorderPad0 :
  case BorderPad : 
  case BorderMirror :
    *pTimeIdMin     = timeIdMin-shiftMax;
    *pTimeIdMax     = timeIdMax-shiftMin;
    break;
  case BorderPeriodic :
    // Case when all changes are in the first 'copy' of the signal
    if(timeIdMax <= stft->signalSize-1) {
      *pTimeIdMin      = timeIdMin-shiftMax;
      *pTimeIdMax      = timeIdMax-shiftMin;
      break;
    }
    // Otherwise : the first part that need to be updated
    *pFlagTwoLoops = YES;
    *pTimeIdMin    = 0-shiftMax;
    *pTimeIdMax    = (timeIdMax%stft->signalSize)-shiftMin;
    *pTimeIdMin1   = timeIdMin-shiftMax;
    *pTimeIdMax1   = stft->signalSize-1-shiftMin;
    break;
  default :
      Errorf("ComputeUpdateLimits : border type %s not treated yet",BorderType2Name(stft->borderType));
  }

  /* We consider only the locations on the  time-freq grid */
  QuantizeRangeLarge(*pTimeIdMin,*pTimeIdMax,stft->tRate,
		     pTimeIdMin,pTimeIdMax);	
  *pTimeIdMin     = MAX(0,*pTimeIdMin);
  *pTimeIdMax     = MIN((stft->nFrames*stft->tRate)-stft->tRate,*pTimeIdMax);

  if(*pFlagTwoLoops) {
    QuantizeRangeLarge(*pTimeIdMin1,*pTimeIdMax1,stft->tRate,
		       pTimeIdMin1,pTimeIdMax1);	
    *pTimeIdMin1     = MAX(0,*pTimeIdMin1);
    *pTimeIdMax1     = MIN((stft->nFrames*stft->tRate)-stft->tRate,*pTimeIdMax1);
  }
}


void AddStft(STFT stftResult,STFT stftAdded,unsigned long timeIdMin,unsigned long timeIdMax)
{
  int timeIdMinStft,timeIdMaxStft;
  int timeIdMinStft2,timeIdMaxStft2;
  char flagTwoLoops;

  unsigned long timeId;
  unsigned long f;
  LWFLOAT *real1,*real2,*imag1=NULL,*imag2;

  CheckStftNotEmpty(stftResult);
  CheckStftNotEmpty(stftAdded);
  CheckSameStftStruct(stftResult,stftAdded);
  ComputeUpdateLimits(stftResult,timeIdMin,timeIdMax,
		      &timeIdMinStft,&timeIdMaxStft,
		      &flagTwoLoops,&timeIdMinStft2,&timeIdMaxStft2);
  // First time Loop 
  for (timeId = timeIdMinStft; 
       timeId <= timeIdMaxStft; 
       timeId += stftResult->tRate) {
    if(stftResult->type==ComplexStft) {
      GetStftData(stftResult,timeId,&real1,&imag1);
      GetStftData(stftAdded,timeId,&real2,&imag2);
    } else {
      GetStftData(stftResult,timeId,&real1,NULL);
      GetStftData(stftAdded,timeId,&real2,NULL);
    }
    for(f = 0; f < stftResult->nSubBands; f++) real1[f]+=real2[f];
    if(imag1) {
      for(f = 0; f < stftResult->nSubBands; f++) imag1[f]+=imag2[f];
    }
  }
  // Second loop if needed
  if(flagTwoLoops) Errorf("????");
} 

void ZeroStft(STFT stft,unsigned long timeIdMin,unsigned long timeIdMax)
{
  int timeIdMinStft,timeIdMaxStft;
  int timeIdMinStft2,timeIdMaxStft2;
  char flagTwoLoops;

  unsigned long timeId;
  unsigned long f;
  LWFLOAT *real1,*imag1=NULL;

  CheckStftNotEmpty(stft);
  ComputeUpdateLimits(stft,timeIdMin,timeIdMax,
		      &timeIdMinStft,&timeIdMaxStft,
		      &flagTwoLoops,&timeIdMinStft2,&timeIdMaxStft2);
  // First time Loop 
  for (timeId = timeIdMinStft; timeId <= timeIdMaxStft; timeId += stft->tRate) {
    if(stft->type==ComplexStft) {
      GetStftData(stft,timeId,&real1,&imag1);
    } else {
      GetStftData(stft,timeId,&real1,NULL);
    }
    for(f = 0; f < stft->nSubBands; f++) real1[f]=0.0;
    if(imag1!=NULL) {
      for(f = 0; f < stft->nSubBands; f++) imag1[f]=0.0;
    }
  }
  // Second loop if needed
  if(flagTwoLoops) Errorf("????");
}