atom.c

LastWave

/*..........................................................................*/
/*                                                                          */
/*      L a s t W a v e    P a c k a g e 'stft' 2.1                         */
/*                                                                          */
/*      Copyright (C) 2000 Remi Gribonval, Emmanuel Bacry and Javier Abadia */
/*      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             */
/*                                                                          */
/*..........................................................................*/


/***********************************************************************************/
/*
 * 	Allocation and Desallocations of atom, management of the '&atom' variable.
 */
/************************************************************************************/

#include "lastwave.h"
#include "int_fsilist.h"
#include "atom.h"

/**********************************/
/*
 * 	ATOM VARIABLES
 */
/**********************************/

char *atomType = "&atom";
/*
 * Answers to the different print messages
 */
 
void ShortPrintAtom(ATOM atom)
{
  PrintAtom(atom,YES);
}

char *ToStrAtom(ATOM atom, char flagShort)
{
  static char str[30];
  sprintf(str,"<&atom;%p>",atom);
  return(str);
}

void PrintInfoAtom(ATOM atom)
{
  PrintAtom(atom,NO);
}

ATOM TNewAtom(void)
{
  ATOM atom;
  
  atom = NewAtom();
  TempValue(atom);
  return(atom);
}

/*******************************************************************/
/*
 * 	Conversion between "sample" coordinates
 *	and "physical" coordinates  using sample-frequency and GABOR_MAX_FREQID
 */
/*******************************************************************/
LWFLOAT LW_TFConvertUnit(const void *content,LWFLOAT value,char unitType,char convertType)
{
  ATFCONTENT tfContent = (ATFCONTENT) content;
  
  // Some checkings
  if(content == NULL)
    Errorf("LW_TFConvertUnit : NULL content");
  if(convertType!= LW_TF_Id2RealConvert && convertType != LW_TF_Real2IdConvert)
    Errorf("LW_TFConvertUnit : (Weird) unknown conversion type %d",convertType);
  if(tfContent->dx <= 0.0)
    Errorf("LW_TFConvertUnit : bad dx %g",tfContent->dx);
  // Conversion
  switch(unitType) {
    // Time
  case LW_TF_TimeUnit :
    if(convertType == LW_TF_Id2RealConvert) 
      return(value*tfContent->dx+tfContent->x0);
    else
      return((value-tfContent->x0)/tfContent->dx);
    // Frequency
  case LW_TF_FreqUnit :
    if(convertType == LW_TF_Id2RealConvert) 
      return(value/(tfContent->dx*GABOR_MAX_FREQID));
    else
      return(value*GABOR_MAX_FREQID*tfContent->dx);
    // Chirp
  case LW_TF_ChirpUnit :
    if(convertType == LW_TF_Id2RealConvert) 
      return(value/(tfContent->dx*tfContent->dx*GABOR_MAX_FREQID));
    else
      return(value*GABOR_MAX_FREQID*tfContent->dx*tfContent->dx);
  default : 
    Errorf("LW_TF_ConvertUnit : (Weird) unknown unit type %d",unitType);
  }
  // Should never be reached but ensures that the function returns a value
  return(0);
}

/**************************************/
/* 
 *	ALLOCATIONS
 */
/**************************************/

/***********************************/
/*
 *           Checkings 
 */
/***********************************/

void CheckAtom(const ATOM atom)
{
  CheckTFContent(atom);
  // Window 
  if(!INRANGE(STFT_MIN_WINDOWSIZE,atom->windowSize,STFT_MAX_WINDOWSIZE))
    Errorf("CheckAtom : bad windowSize %d",atom->windowSize);
  if(atom->windowSize>atom->signalSize)
    Errorf("CheckAtom : bad windowSize %d compared to signalSize %d",atom->windowSize,atom->signalSize);
  if(!WindowShapeIsOK(atom->windowShape))
    Errorf("CheckAtom : bad windowShape %d",atom->windowShape);
  // Time location
  if(atom->timeId != (int) atom->timeId)  Errorf("CheckAtom : timeId is not int (not implemented yet)");      
  if(!INRANGE(0,atom->timeId,atom->signalSize-1))
    Errorf("CheckAtom : bad timeId %g not in [0 %d]",atom->timeId,atom->signalSize-1);
  // Freq location : TODO remove this restriction to allow pitch shifting ?
  // Or restrict to GABOR_NYQUIST_FREQID after changing the "AutoInnerProd" function ?
  if(!INRANGE(0,atom->freqId,GABOR_MAX_FREQID))
    Errorf("CheckAtom : bad freqId %g not in [0 %d]",atom->freqId,GABOR_MAX_FREQID);
  // Chirp `location'
  if(!INRANGE(-GABOR_MAX_CHIRPID,atom->chirpId,GABOR_MAX_CHIRPID))    
    Errorf("CheckAtom : bad |chirpId|=|%g| > %d",atom->chirpId,GABOR_MAX_CHIRPID);

  // Case of DC/Nyquist atoms
  if(atom->freqId == 0.0 || atom->freqId == GABOR_NYQUIST_FREQID) {
    if(atom->chirpId != 0.0)
      Errorf("CheckAtom : DC/Nyquist (freqId %g) is chirped at %g",atom->freqId,atom->chirpId);
    if(atom->coeffI != 0.0)
      Errorf("CheckAtom : DC/Nyquist (freqId %g) has coeffI %g",atom->freqId,atom->coeffI);
  }

  if(atom->coeff2 < 0) Errorf("CheckAtom : coeff2 %g < 0!",atom->coeff2);
  if(atom->realGG*atom->realGG+atom->imagGG*atom->imagGG>1)
    Errorf("CheckAtom : bad GG (%g,%g)",atom->realGG,atom->imagGG);
}

void CheckAtomReal(const ATOM atom)
{
  CheckAtom(atom);
  if(atom->flagGGIsSet == NO)
    Errorf("CheckAtomReal : atom GG is not set");
}

/*  Allocation/ Initialization of ATOM */
static void InitAtom(ATOM atom)
{
  InitTFContent(atom);

  // Window shape
  atom->windowShape = STFT_DEFAULT_WINDOWSHAPE;
  atom->windowSize  = STFT_DEFAULT_WINDOWSIZE;  
  // Time-freq location
  atom->timeId  = 0.0;
  atom->freqId  = 0.0;
  atom->chirpId = 0.0;
  
  // This is the right value of gg when freqId==0
  atom->realGG   = 1.0;
  atom->imagGG   = 0.0;
  atom->flagGGIsSet = NO;

  // The complex coeff is zero
  atom->coeffR      = 0.0;
  atom->coeffI      = 0.0;
  // The squared amplitude and phase 
  atom->coeff2   = 0.0;
  atom->cosPhase = 1.0;
  atom->sinPhase = 0.0;

}

ATOM NewAtom(void)
{
  ATOM atom;
  
#ifdef DEBUGALLOC
  DebugType = "Atom";
#endif
  
  atom = (ATOM) Malloc(sizeof(Atom));
  InitValue(atom,&tsAtom); 
  InitAtom(atom);
  return(atom);
}

void ClearAtom(ATOM atom)
{
  InitAtom(atom);
}

ATOM DeleteAtom(ATOM atom)
{
  if(atom == NULL) Errorf("DeleteAtom : NULL atom");
  
  if (atom->nRef==0) Errorf("DeleteAtom : Weird Error : trying to delete a temporary atom\n");
  
  RemoveRefValue(atom);
  if (atom->nRef > 0) return(NULL);
  
#ifdef DEBUGALLOC
  DebugType = "Atom";
#endif
  
  Free(atom);
  return(NULL);
}

// Makes a copy of all fields of an atom
ATOM CopyAtom(const ATOM in,ATOM out)
{
  if(in == NULL) return(NULL);
  if(out == NULL) out = NewAtom();
  if(in == out) return(out);

  CopyFieldsTFContent(in,out);
  out->windowShape = in->windowShape;
  out->windowSize  = in->windowSize;  

  out->timeId   = in->timeId;
  out->freqId   = in->freqId;
  out->chirpId  = in->chirpId;
  
  out->realGG     = in->realGG;
  out->imagGG     = in->imagGG;  
  out->flagGGIsSet  = in->flagGGIsSet;

  out->coeffR = in->coeffR;
  out->coeffI = in->coeffI;
  out->coeff2   = in->coeff2;
  out->cosPhase = in->cosPhase;
  out->sinPhase = in->sinPhase;

  return(out);
}

/* Printing utility */
void PrintAtom(const ATOM atom,char flagShort)
{
  CheckAtom(atom);
  if(flagShort) {
    Printf("(s,t,f,c) = (%d,%g,%g,%g), ",
	   atom->windowSize,TimeId2Time(atom,atom->timeId),FreqId2Freq(atom,atom->freqId),ChirpId2Chirp(atom,atom->chirpId));
    if(atom->flagGGIsSet)
      Printf("coeff2=%g phase=(%g,%g)\n",atom->coeff2,atom->cosPhase,atom->sinPhase);
    else
      Printf("coeff (%g,%g)",atom->coeffR,atom->coeffI);
  }
  else {
    Printf("Scale  : %g (%d)\n",atom->dx*atom->windowSize,atom->windowSize,atom->dx);
    Printf("Time   : %g (%g)\n",TimeId2Time(atom,atom->timeId),atom->timeId);
    Printf("Freq   : %g (%g)\n",FreqId2Freq(atom,atom->freqId),atom->freqId);
    Printf("Chirp  : %g (%g)\n",ChirpId2Chirp(atom,atom->chirpId),atom->chirpId);
    if(atom->flagGGIsSet) {
#ifdef ATOM_ADVANCED
    Printf("  correlation with complex conjugate  (%g,%g)\n",atom->realGG,atom->imagGG);
#endif // ATOM_ADVANCED
      Printf("Coeff2 : %g\n",atom->coeff2);
      Printf("Phase  : (%g,%g)\n",atom->cosPhase,atom->sinPhase);
    } 
    Printf("Coeff  : (%g,%g)\n",atom->coeffR,atom->coeffI);
  }
}

/* returns YES iff the time-support of the atoms intersect
 * and [*pTimeIdMin,*pTimeIdMax] is the intersection of the supports
 */
char AtomsIntersect(const ATOM atom1,const ATOM atom2,long *pTimeIdMin,long *pTimeIdMax)
{
  int timeIdMin1,timeIdMin2,timeIdMax1,timeIdMax2;
  long timeIdMin,timeIdMax;
  
  /* Checking arguments */
  CheckAtom(atom1);
  CheckAtom(atom2);
  
  /* Computing supports and their intersection */
  ComputeWindowSupport(atom1->windowSize,atom1->windowShape,atom1->timeId,&timeIdMin1,&timeIdMax1);
  ComputeWindowSupport(atom2->windowSize,atom2->windowShape,atom2->timeId,&timeIdMin2,&timeIdMax2);
  
  timeIdMin = MAX(timeIdMin1,timeIdMin2);
  timeIdMax = MIN(timeIdMax1,timeIdMax2);
  
  if(pTimeIdMin != NULL)
    *pTimeIdMin = timeIdMin;
  if(pTimeIdMax != NULL)
    *pTimeIdMax = timeIdMax;
  
  return(timeIdMin<=timeIdMax);
}

/*
 *
 * Computes (realGG,imagGG) : The real and imaginary parts
 * of the integral of the atom with its conjugate
 *
 * We guarantee that if the atom has symmetries one has imagGG = 0
 * Namely : if the chirp is zero and the window is symmetric (e.g., Gaussian)
 *
 */
static void SetAtomGG(ATOM atom)
{
  static ATOM copy = NULL;
  
  /* The signals to build the atom in */
  static SIGNAL  atomSignalR = NULL;
  static SIGNAL  atomSignalI = NULL;
  
  LWFLOAT realGG,imagGG;
  
  /* Utility */
  char flagAsym = NO;
  /* Checkings */
  CheckAtom(atom);
  // Case where the GG is tabulated
  if(atom->chirpId == 0.0 && atom->freqId == (int) atom->freqId) {
    if(ComputeWindowGG(atom->windowShape,atom->windowSize,atom->freqId,&realGG,&imagGG)) {
      atom->realGG = realGG;
      atom->imagGG = imagGG;
      atom->flagGGIsSet = YES;
      return;
    }
  }
  // Else, do a numeric or a fast analytic computation 
  AutoAtomInnerProduct(atom,NO,&realGG,&imagGG);
  atom->realGG = realGG;
  atom->imagGG = imagGG;
  atom->flagGGIsSet = YES;
}

void TouchAtomGG(ATOM atom)
{
  atom->flagGGIsSet = NO;
}



// TODO : explain 
void ComputeRealPhaseEnergy(LWFLOAT coeffR,LWFLOAT coeffI,LWFLOAT realGG,LWFLOAT imagGG,
			    LWFLOAT *pPhase,LWFLOAT *pCosPhase,LWFLOAT *pSinPhase,
			    LWFLOAT *pEnergy)
{
  LWFLOAT innerProd2;
  LWFLOAT real,imag;
  LWFLOAT norm;

  // Basic checking
  if(realGG*realGG+imagGG*imagGG>1) Errorf("ComputeRealPhaseEnergy : (Weired) bad GG (%g %g)",realGG,imagGG);
  if(pPhase==NULL && (pCosPhase!=NULL || pSinPhase!=NULL))
    Errorf("ComputeRealPhaseEnergy : bad phase pointers");
  if(pPhase!=NULL && (pCosPhase==NULL || pSinPhase==NULL))
    Errorf("ComputeRealPhaseEnergy : bad phase pointers");

  // The squared magnitude of the complex innerproduct
  innerProd2 	= coeffR*coeffR+coeffI*coeffI;

  // The simplest case : when the result is zero!
  if (innerProd2 == 0)  {
    if(pEnergy) 
      *pEnergy 	= 0.0;
    if(pPhase) {
      *pCosPhase= 1.0;
      *pSinPhase= 0.0;
      *pPhase 	= 0.0;
    }
    return;
  }
  
  // The special case where (realGG,imagGG)==(1,0)
  // corresponds to freqId = 0,Nyquist : 
  // -> no phase, just a sign
  if(realGG==1.0) {
    // Checking :
    if(imagGG != 0.0 || coeffI != 0.0) Errorf("ComputeRealPhaseEnergy : non zero imaginary part when realGG==1.0!");

    // Setting energy
    if(pEnergy)
      *pEnergy = innerProd2;
    // Setting phase
    if(pPhase) {