stft_tabulate.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 and J.McKelvie         */
/*      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             */
/*                                                                          */
/*..........................................................................*/



#include "lastwave.h"

#include "stft.h"


/*
 * To get rid of
 */
 
 
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

#define FD LWFLOAT

static void FourierTransform(FD *data,int nn,int fftSign)
{
	int n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	LWFLOAT tempr,tempi;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(data[j],data[i]);
			SWAP(data[j+1],data[i+1]);
		}
		m=n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=2*mmax;
		theta=6.28318530717959/(fftSign*mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {
			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*data[j]-wi*data[j+1];
				tempi=wr*data[j+1]+wi*data[j];
				data[j]=data[i]-tempr;
				data[j+1]=data[i+1]-tempi;
				data[i] += tempr;
				data[i+1] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}

#undef SWAP

/* Get the real part of the fftData and put it in the signal 'out' */

static void GetRealPart(FD *fftData,int size,SIGNAL out)
{
   int i;

   SizeSignal(out,size,YSIG);
   for (i=0;i<size;i++)
	   out->Y[i] = fftData[2*i+1];
}


/* Get the imaginary part of the fftData and put it in the signal 'out' */

static void GetImagPart(FD *fftData,int size,SIGNAL out)
{
   int i;
   
   SizeSignal(out,size,YSIG);
   for (i=0;i<size;i++)
	  out->Y[i] = fftData[2*i+2];
}

/***********************************************************************/
/* Make FFT-form signal from realPart and imagPart                     */
/* Add zeros if the size != newSize (newSize must be a power of 2      */
/* If imagPart == NULL then it is considered to be made of zero        */
/***********************************************************************/

static void MakeFFTData(FD *fftData,LWFLOAT *realPart,LWFLOAT *imagPart,int size,int newSize)
{
   int   i;

   if (imagPart != NULL) {     
   	for (i=0;i<newSize;i++)
       	 if (i>=size)
          	fftData[2*i+1] = fftData[2*i+2] = 0.0;
         else {
           	fftData[2*i+1] = realPart[i];
           	fftData[2*i+2] = imagPart[i];
         }
     }
    else {
      for (i=0;i<newSize;i++)
       	 if (i>=size)
          	fftData[2*i+1] = fftData[2*i+2] = 0.0;
         else {
           	fftData[2*i+1] = realPart[i];
           	fftData[2*i+2] = 0.0;
         }
     }  
}

/***********************************************************/
/*                                                         */
/*         Compute the fft transform of a SIGNAL           */
/*             (inImag can be NULL)                        */
/*                                                         */
/***********************************************************/


/* fft :  SIGNAL --> float *          */
/* fftOut must be of size 2*newSize+2 */
static void Fft1_(SIGNAL inReal,SIGNAL inImag,FD *fftOut, int newSize, int fftSign)
{

  if (inImag == NULL) 
   	  MakeFFTData(fftOut,inReal->Y,NULL,inReal->size,newSize);
  else
      MakeFFTData(fftOut,inReal->Y,inImag->Y,inReal->size,newSize);
    	  
  FourierTransform(fftOut,newSize,fftSign);
}
   

/* fft :  SIGNAL --> SIGNAL          */
static void Fft1(SIGNAL inReal,SIGNAL inImag,SIGNAL outReal,SIGNAL outImag,int fftSign,char flagShift)
{
   int size2, newSize,size;
   FD *fftData,val;
   int i;
   
   size2 = (int) (0.5+(log((double) inReal->size)/log(2.)));
   newSize = 1 << size2;
   
   /* allocation of the result */
   fftData = (FD *) Malloc(sizeof(FD)*(2*newSize+2));
   
   Fft1_(inReal,inImag,fftData,newSize,fftSign);  
   
   GetRealPart(fftData,newSize,outReal);
   GetImagPart(fftData,newSize,outImag);
   
   if (fftSign == -1) 
   	for(i=0;i<outReal->size;i++) {
   		outReal->Y[i] /= newSize;
   		outImag->Y[i] /= newSize;
   	}
   
   if (!flagShift) {
     outImag->dx = outReal->dx = 2*M_PI/outImag->size; 
     outImag->x0 = outReal->x0 = 0;
   } 
   else {
     size = outImag->size;
     for (i=0;i<size/2;i++ ) {
       val = outImag->Y[i];
       outImag->Y[i] = outImag->Y[size/2+i];
       outImag->Y[size/2+i] = val;
       val = outReal->Y[i];
       outReal->Y[i] = outReal->Y[size/2+i];
       outReal->Y[size/2+i] = val;
     }
     outImag->dx = outReal->dx = 1/(outImag->size*inReal->dx); 
     outImag->x0 = outReal->x0 = -1/(2*inReal->dx);
   }
   
   
   Free(fftData);
}


//
// Tabulation of some data for the Stft package
//

//
// THE STATIC VARIABLES BELOW ARE INITIALIZED IN InitStftTabulation 
// AND FREED IN FreeStftTabulation
//

// The number of values of 'windowSize' for which some tabulation occurs
// and an array that maps i->windowSize for 0<= i < nTabWindowSizes
static unsigned short nTabWindowSizes = 0;
static unsigned long* stftTabWindowSizes = NULL;
// The function to convert from a value to an index 0<=i<nTabWindowSizes
// If 'windowSize' is not tabulated we return nTabWindowSizes;
static unsigned short GetIndexWindowSize(unsigned long windowSize) {
  unsigned short i;
  for(i=0;i<nTabWindowSizes;i++) {
    if(stftTabWindowSizes[i]==windowSize) return(i);
  }
  return(i);
}

// Memory is allocated for these signals in InitStftTabulations
static SIGNAL** stftTabGGI;
static SIGNAL** stftTabGGR;

/**********************************************/
/*   TODO : explain
 *   Tabulation of the GG
 *   (realGG[o-oMin]->Y[freqId],imagGG[o-oMin]->Y[freqId]) 
 *   is the inner product <g,\overline{g}> for the atom g
 *   at octave 'o' and frequency freqId.
 *
 *   -For symmetric 'windowShape', 'imagGG' is always zero so 
 *   *pImagGG is set to NULL
 *
 *   -(realGG,imagGG) at frequency freqId is the ?? as at freqId
 *    (GABOR_NYQUIST_FREQID-freqId) so we don't need to store the frequencies for
 *    freqId > FreqNyquist/2.
 *
 *   -the values |(realGG,imagGG)| <= CCATOM_PRECISION are not stored.
 *
 * -----------------------------------------------------------------
 *   -at the zero frequency
 *	realGG[o-oMin]->Y[0] == 1.0
 *	imagGG[o-oMin]->Y[0] == 0.0
 *
 ***********************************************************/

static void StftTabulateGG(char windowShape)
{ /* Declaration of variables used in 'for' loops below */
  unsigned short i;
  unsigned long windowSize;
  SIGNAL window,paddedWindow,shiftedWindow;
  char flagAsymWindow = NO;
  // TODO : explain
  unsigned long maxWindowShape;

  SIGNAL *tab1GGR  = NULL;
  SIGNAL *tab1GGI  = NULL;
  SIGNAL GGR      = NULL;
  SIGNAL GGI      = NULL;

  LWFLOAT realGG,imagGG;
  unsigned long freqId,n;

  // TODO : remove !
  LWFLOAT precision = CCATOM_PRECISION*CCATOM_PRECISION;

  /* Checking arguments */
  if(!WindowShapeIsOK(windowShape)) Errorf("StftTabulateGG : bad windowShape %d",windowShape);
  // Quit without dong anything if already tabulated
  if(stftTabGGR[windowShape]!=NULL) return;

  flagAsymWindow = GetFlagAsymWindowShape(windowShape);

  // TODO : try to load from file here

  /* Allocating */
  tab1GGR = stftTabGGR[windowShape] = (SIGNAL*) Calloc(nTabWindowSizes,sizeof(SIGNAL));
  tab1GGI = stftTabGGI[windowShape] = (SIGNAL*) Calloc(nTabWindowSizes,sizeof(SIGNAL));
  paddedWindow  = NewSignal();
  shiftedWindow = NewSignal();
  SizeSignal(paddedWindow,GABOR_MAX_FREQID,YSIG);
  SizeSignal(shiftedWindow,GABOR_MAX_FREQID,YSIG);
  
  //DEBUG
  Printf("Tabulating GG for '%s' windows\n",WindowShape2Name(windowShape));Flush();

  // Allocating the signals and filling them
  for(i=0; i<nTabWindowSizes; i++) {
    windowSize = stftTabWindowSizes[i];
    GetTabWindow(windowShape,windowSize,&window);
    // TODO : CLARIFY AND EXPLAIN THIS 
    maxWindowShape = GetMaxWindowShape(windowShape,windowSize);
    //DEBUG
//    Printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
    Printf("***tabulating size %6d\n",windowSize);Flush();
    //DEBUG
    Printf("--maxWindowShape=%6d\n",maxWindowShape);Flush();
    
    // Initializing squared window
    for(n=0; n<window->size; n++)    paddedWindow->Y[n] = window->Y[n]*window->Y[n];
    // Zero padding
    for(; n<paddedWindow->size; n++) paddedWindow->Y[n] = 0.0;
    // Rotation to take into account the maxWindowShape : 
    // the window must be 'centered' at maxWindowShape 
    for(n = 0; n<shiftedWindow->size; n++) 
      shiftedWindow->Y[n] = paddedWindow->Y[(n+maxWindowShape)%paddedWindow->size];		

    //
    // Allocating and tabulating the resultin GG signals
    //