stft_complex.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 COMPLEX Short Time Fourier Transform :
 *
 *		COMPUTATION and UPDATE
 */
/****************************************************/
#include "lastwave.h"

#include "stft.h"


/*
 * FFT : To get rid of (use a library!)
 */
 
 
#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);
}


/************************************************************************
 *
 *    Multiply signal with the window and FFT
 *   (Sum_n  signal(n+shift) * window(n) * exp(2*i*PI*k*fRate*n/GABOR_MAX_FREQID)
 *    
 *
 ************************************************************************/
void FilterMultiplyFft(SIGNAL signal,char borderType,SIGNAL window,char windowShape,unsigned long timeId, unsigned long fRate,
		       SIGNAL resultR,SIGNAL resultI)
{
  unsigned long maxWindowShape;
  unsigned long paddFftSize;
  unsigned long i;
  static SIGNAL timeSignal = NULL;
  
  /* Checking arguments */
  if(signal == NULL) Errorf("FilterMultiplyFft : NULL signal");
  if(window == NULL) Errorf("FilterMultiplyFft : NULL window");
  if(window->size%2 != 0) Errorf("FilterMultiplyFft   : window size %d should be even",window->size);
  if(GABOR_MAX_FREQID%fRate != 0) Errorf("FilterMultiplyFft   : bad fRate %d does not divide GABOR_MAX_FREQID",fRate,GABOR_MAX_FREQID);
  
  /* Allocation if necessary */
  if (timeSignal == NULL) timeSignal = NewSignal();
  
  /* 
   * Compute the size of the (possibly zero padded) 
   * signal I want to perform the FFT on :
   * I need at least 'GABOR_MAX_FREQID/fRate' bins.
   */
  paddFftSize = MAX(window->size,GABOR_MAX_FREQID/fRate);
  
  /* 
   * Extract just the piece of the signal that we need,
   * with treatment of the border effects.
   * 
   * Checked correct on the 21/02/2001 by R.Gribonval
   *
   * timeId-maxWindowShape is the index of the first 
   * point of the window (which value is set to zero)
   * This is also the first point of the signal to extract.
   */
  maxWindowShape = GetMaxWindowShape(windowShape,window->size); 
  ExtractSig(signal,timeSignal,borderType,timeId-maxWindowShape,paddFftSize);
  
  /* Multiplication by window */
  for(i=0; i < window->size; i++) timeSignal->Y[i] *= window->Y[i];
    
  /* Zero padding if necessary */
  for(; i < paddFftSize; i++) timeSignal->Y[i]  = 0.0;
  
  /* Rotation to take into account the maxWindowShape : 
   * the window must be 'centered' at maxWindowShape */
  SizeSignal(resultR,paddFftSize,YSIG);
  ZeroSig(resultR);
  for(i=0; i<resultR->size; i++) 
    resultR->Y[i] = timeSignal->Y[(i+maxWindowShape)%timeSignal->size];		
  
  /* Fft */
  SizeSignal(resultI,paddFftSize,YSIG);
  ZeroSig(resultI);
  Fft1(resultR,NULL,resultR,resultI,1,NO);
  
  /* Taking the conjugate value : Why ???? */
  for(i=0; i<resultR->size; i++) resultI->Y[i] *= -1;
  
  /* Properly setting the 'dx' field */
  resultR->dx = GABOR_MAX_FREQID/paddFftSize;
  resultI->dx = GABOR_MAX_FREQID/paddFftSize;
}

void C_FilterMultiplyFft(char **argv)
{
  SIGNAL signal;
  char borderType = STFT_DEFAULT_BORDERTYPE;
  unsigned long fRate;
  char windowShape = STFT_DEFAULT_WINDOWSHAPE;

  SIGNAL window ;
  unsigned long timeId;
  SIGNAL resultR,resultI;
  char opt;
  char *name;
  argv = ParseArgv(argv,tSIGNALI,&signal,tSIGNALI,&window,tINT,&timeId,tINT,&fRate,tSIGNAL,&resultR,tSIGNAL,&resultI,-1);
  
  while((opt = ParseOption(&argv))) {
    switch(opt) {
    case 'b' :
      argv = ParseArgv(argv,tSTR,&name,-1);
      borderType = Name2BorderType(name);
      break;
    default :
      ErrorOption(opt);
    }
  }
  NoMoreArgs(argv);
  FilterMultiplyFft(signal,borderType,window,windowShape,timeId,fRate,resultR,resultI);
}
#ifdef TODO
void C_OverlapAdd(char **argv)
{