stft_complex.c

LastWave

  STFT stft;
  SIGNAL signalR,signalI;
  int timeId;
  LWFLOAT *real,*imag;
  argv = ParseArgv(argv,tSTFT_,NULL,&stft,tSIGNAL,&signalR,tSIGNAL,&signalI,0);
  if(stft == NULL)
    stft = GetStftCur();
  CheckStftComplex(stft);
  SizeSignal(signalR,stft->signalSize,YSIG);
  SizeSignal(signalI,stft->signalSize,YSIG);
  
  for(timeId = 0; 
      timeId < stft->nFrames*stft->tRate;
      timeId += stft->tRate) {
    GetStftData(stft,timeId,&real,&imag);

    MakeFFTData(fftOut,real,imag,stft->nSubBands,newSize);
    FourierTransform(fftOut,newSize,0);
  }
}
#endif // TODO

/*******************************************************************
*     Convolution of a signal with a complex filter
************************ WARNING************************************
* We want to do the convolution with the window reversed
* so that we obtain scalar product (modulo the conjugate)
* Beware of the asymmetric windows !
************************ WARNING*************************************/
void FilterConvol(SIGNAL signal,int borderType,SIGNAL window,int windowShape,int freqId,int tRate,
		  SIGNAL resultR,SIGNAL resultI)
{
  static SIGNAL filterR = NULL;
  static SIGNAL filterI = NULL;
  SIGNAL expikR = NULL;
  SIGNAL expikI = NULL;
  int maxWindowShape;
  int index,ishift;
  int i;
  
  /* Checking arguments */
  if(signal == NULL) Errorf("FilterConvol : NULL signal");
  if(window == NULL) Errorf("FilterConvol : NULL window");
  if(window->size%2 != 0) Errorf("FilterConvol   : window size %d should be even",window->size);

  /* 
   * Build the complex filter
   */
  /* Allocating (only once) */
  if(filterR == NULL) {
    filterR = NewSignal();
    filterI = NewSignal();
  }
  /* Their window =  reversed time window */
  SizeSignal(filterR,window->size,YSIG);
  SizeSignal(filterI,window->size,YSIG);
  for(i=0;i<window->size;i++) {
    index = window->size-1-i;
    filterR->Y[i] = window->Y[index];
    filterI->Y[i] = window->Y[index];
  }
  /*
   * Their modulation = exp(2*i*pi*freqId*(i-ishift)/GABOR_MAX_FREQID)
   * with i-ishift=0 when i is at the maximum of the window
   * of the filter
   * i.e. when i = window->size-1-maxWindowShape
   */
  maxWindowShape = GetMaxWindowShape(windowShape,window->size);
  ishift        = window->size-1-maxWindowShape;

  GetTabExponentials(&expikR,&expikI);
  for(i=0;i<window->size;i++) {
    index = (freqId*(i-ishift+GABOR_MAX_FREQID))%GABOR_MAX_FREQID;
    if(!INRANGE(0,i,expikR->size-1))
      Errorf("FilterConvol : (Weird) internal expikR");
    filterR->Y[i] *= expikR->Y[index];
    filterI->Y[i] *= expikI->Y[index];
  }

  /* The maximum of the filter is at ishift. 
   * We shift it by -ishift to put the maximum at zero */
  filterR->x0 = -ishift;
  filterI->x0 = -ishift;
  
  ConvSig(signal,filterR,resultR,borderType,CV_UNDEFINED);
  ConvSig(signal,filterI,resultI,borderType,CV_UNDEFINED);
}

void C_FilterConvol(char **argv)
{
  SIGNAL signal;
  char borderType = STFT_DEFAULT_BORDERTYPE;
  int octave;
  static SIGNAL window = NULL;
  int windowShape  = STFT_DEFAULT_WINDOWSHAPE;
  int freqId,tRate;

  SIGNAL resultR,resultI;
  
  char opt;
  char *name;
  
  argv = ParseArgv(argv,tSIGNALI,&signal,tINT,&octave,tINT,&freqId,tINT,&tRate,
		   tSIGNAL,&resultR,tSIGNAL,&resultI,-1);
  
  while((opt = ParseOption(&argv))) {
    switch(opt) {
    case 'b' :
      argv = ParseArgv(argv,tSTR,&name,-1);
      borderType = Name2BorderType(name);
      break;
    case 'w' :
      argv = ParseArgv(argv,tSTR,&name,-1);
      windowShape = Name2WindowShape(name);
      break;
    default :
      ErrorOption(opt);
    }
  }
  NoMoreArgs(argv);

  /* Allocate once, if necessary */
  if(window == NULL) {
    window  = NewSignal();
  }

  GetTabWindow(windowShape,1<<octave,&window);
  FilterConvol(signal,borderType,window,windowShape,freqId,tRate,resultR,resultI);
}

 
/**************************************************************/
/*
 * Update of the COMPLEX stft data, two methods :
 *
 *	'Time' : we loop on time, within a SPECIFIED TIME RANGE
 *		     and update all frequencies with a 
 *			'filter-multiply + FFT' algorithm
 *	'Freq' : we loop on ALL frequencies and update with a 
 *                    CONVOLUTION (Direct or FFT = optimized?)
 */
/**************************************************************/


/*
 * Very helpful function to set the values of a complex stft 
 * at a given time 'timeId', for all frequencies.
 * 'freqSignal' is a signal of 'total size' (size*dx) 'GABOR_MAX_FREQID'
 * representing the values of the complex spectrogram
 * at different frequencies (from 0 to Nyquist=GABOR_MAX_FREQID/2 with a step freqSignal->dx)
 */
static void StftComplexSetFreq(STFT stftComplex,SIGNAL freqSignalR,SIGNAL freqSignalI,int timeId)
{
    int freqId,fC,sigScale;
    LWFLOAT *innerProdR,*innerProdI;

    /* Checking arguments */
    CheckStftComplex(stftComplex);
    if(!INRANGE(0,timeId,stftComplex->signalSize-1))
      Errorf("StftComplexSetFreq() : bad timeId %d not in [0 %d]",timeId,stftComplex->signalSize-1);
    if(freqSignalR == NULL || freqSignalI == NULL)
      Errorf("StftComplexSetFreq() : NULL freqSignalR or freqSignalI");
    if(freqSignalR->size != freqSignalI->size)
      Errorf("StftComplexSetFreq() : freqSignalR and freqSignalI have different size");
    if(freqSignalR->dx != freqSignalI->dx)
      Errorf("StftComplexSetFreq() : freqSignalR and freqSignalI have different dx");
    sigScale = (int) freqSignalR->dx;
    if (sigScale != stftComplex->fRate)
	  Errorf("StftComplexSetFreq() : fRate %d != signal dx %d\n", stftComplex->fRate,sigScale);
    if (freqSignalR->size*freqSignalR->dx != GABOR_MAX_FREQID)
      Errorf("StftComplexSetFreq() : freqSignal total size %g must be GABOR_MAX_FREQID %d\n",
	     freqSignalR->size*freqSignalR->dx,GABOR_MAX_FREQID);

    /* Getting the complex data location */
    GetStftData(stftComplex,timeId,&innerProdR,&innerProdI);
  
    /* We first deal with the frequencies 0 */
    innerProdR[0] 	= freqSignalR->Y[0];
    innerProdI[0] 	= 0.0;	/* The inner-product is real */

    /* Then we deal with frequency Nyquist (if exists) */
    if (stftComplex->nSubBands >= 2) {
      fC = GABOR_NYQUIST_FREQID/stftComplex->fRate;
      // TODO : check it is correct !!!
      innerProdR[fC] 	= freqSignalR->Y[GABOR_NYQUIST_FREQID/sigScale];
      innerProdI[fC] 	= 0.0;  /* The inner-product is real */
    }

    /* Then we deal with the other frequencies */
    for(freqId = stftComplex->fRate, fC = 1; fC < stftComplex->nSubBands-1; freqId += stftComplex->fRate,fC++) {
      innerProdR[fC] 	= freqSignalR->Y[freqId/sigScale];
      innerProdI[fC] 	= freqSignalI->Y[freqId/sigScale];
    }
} 

/**************************************************************
 *
 * Computes the values of the COMPLEX stft of a given signal,
 * given that this signal has changed between the time-indexes
 * [timeIdMin, timeIdMax]. 
 * If the borderType is BorderPad or BorderPad0, 
 *   [timeIdMin,timeIdMax] must be a subinterval of [0,signalSize-1]
 * If the borderType is BorderPeriodic
 *   [timeIdMin,timeIdMax] must be a subinterval of [0,2*signalSize-1]
 **************************************************************/

void UpdateStftComplex(STFT stftComplex,SIGNAL signal,int timeIdMin,int timeIdMax)
{    
    static SIGNAL freqSignalR = NULL;
    static SIGNAL freqSignalI = NULL;

    char flagTwoLoops = NO;
    int timeIdMinStft,timeIdMaxStft;
    int timeIdMinStft2,timeIdMaxStft2;

    int timeId;

    SIGNAL window = NULL;
    
    int octave;

    /* Checking arguments */
    CheckStftComplex(stftComplex);
    if(signal == NULL) Errorf("UpdateStftComplex : NULL signal");
    
    /* Allocation if necessary */
    if(freqSignalR == NULL) {
	  freqSignalR = NewSignal();	
	  freqSignalI = NewSignal();
    }

    ComputeUpdateLimits(stftComplex,timeIdMin,timeIdMax,
			&timeIdMinStft,&timeIdMaxStft,
			&flagTwoLoops,&timeIdMinStft2,&timeIdMaxStft2);
    /**************************
     *
     * Actual Computation
     *
     **************************/
    
    octave = (int) (log(stftComplex->windowSize)/log(2.0)+0.5);
    if(1<<octave != stftComplex->windowSize)
      Errorf("UpdateStftComplex : windowSize is not a power of two!");
    GetTabWindow(stftComplex->windowShape,stftComplex->windowSize,&window);
    // First time loop
    for (timeId = timeIdMinStft;
	 timeId <= timeIdMaxStft; 
	 timeId += stftComplex->tRate) {
      // Gets the FFT of signal(timeId+n)*window(n) in a complex signal.
      FilterMultiplyFft(signal,stftComplex->borderType,window,stftComplex->windowShape,timeId,stftComplex->fRate,freqSignalR,freqSignalI);
      // Set the complex spectrogram data
      StftComplexSetFreq(stftComplex,freqSignalR,freqSignalI,timeId);
    }
    // Second loop if needed
    if(flagTwoLoops) {
      for (timeId = timeIdMinStft2; 
	   timeId <= timeIdMaxStft2; 
	   timeId += stftComplex->tRate) {
	FilterMultiplyFft(signal,stftComplex->borderType,window,stftComplex->windowShape,timeId,stftComplex->fRate,freqSignalR,freqSignalI);
	StftComplexSetFreq(stftComplex,freqSignalR,freqSignalI,timeId);
      }
    }
    
    /* Now the COMPLEX spectrogram is up to date */
    stftComplex->flagUpToDate = YES;
}

/************************************************************************/
/*
 *			 FREQUENCY LOOP 
 *
 *
 * Computes the values of the complex stft of the given signal,
 * for all time and freq.
 *
 */
/*************************************************************************/
/* Frequency Loop NOT YET RE-IMPLEMENTED */

void UpdateStftComplexFreq(STFT stftComplex,SIGNAL signal)
{
    int freqId;
    static SIGNAL timeSignalR = NULL;
    static SIGNAL timeSignalI = NULL;
    SIGNAL window    = NULL;

    int octave;

    /* Checking arguments */
    CheckStftComplex(stftComplex);
    Errorf("UpdateStftComplexFreq : not re-implemented yet");

    /* Initializing (just once) */
    if(timeSignalR == NULL) {
      timeSignalR = NewSignal();
      timeSignalI = NewSignal();
    }

    octave = (int) (log(stftComplex->windowSize)/log(2.0)+0.5);
    if(1<<octave != stftComplex->windowSize)
      Errorf("UpdateStftComplexFreq : windowSize is not a power of two!");
    GetTabWindow(stftComplex->windowShape,stftComplex->windowSize,&window);

    /* Frequency Loop */
    for (freqId = 0; freqId < stftComplex->fRate*stftComplex->nSubBands; freqId += stftComplex->fRate) {
      FilterConvol(signal,stftComplex->borderType,window,stftComplex->windowShape,
		   freqId,stftComplex->tRate,timeSignalR,timeSignalI);
      Errorf("UpdateStftComplexFreq : subsampling not re-implemented yet");
/*	StftComplexSetTimeData(stftComplex,freqId,timeSignalR,timeSignalI);*/
    }
    /* Now the complex part is up-to-date */
    stftComplex->flagUpToDate = YES;
}

/* EOF */