atom_innerprod.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             */
/*                                                                          */
/*..........................................................................*/


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

/***********************************************/
/* Inner-products between time-frequency atoms */
/***********************************************/

// External functions that perform approximate analytic computation of the inner product
// between two complex atoms for some windowShapes
extern void CCAtomAnalyticGauss(const ATOM atom1,const ATOM atom2,LWFLOAT *pReal,LWFLOAT *pImag);
extern void CCAtomAnalyticFoF(const ATOM atom1,const ATOM atom2,LWFLOAT *pReal,LWFLOAT *pImag);

// The EXACT inner-product between two complex atoms, with numerical computations (it involves building the atoms).
//  The values of coeffR/coeffI/coeff2/realGG/imagGG/cosPhase/sinPhase 
//  do not influence the result.
// TODO : optimize it (it is NOT AT ALL OPTIMIZED!!!)
static void CCAtomNumericInnerProduct(const ATOM atom1,const ATOM atom2,LWFLOAT *pReal,LWFLOAT *pImag)
{
  static SIGNAL atomR1 = NULL;
  static SIGNAL atomI1 = NULL;
  static SIGNAL atomR2 = NULL;
  static SIGNAL atomI2 = NULL;

  // TODO : other border types ?
  char borderType = STFT_DEFAULT_BORDERTYPE;
  long freqId1,freqId2;

  long timeIdMin,timeIdMax;
  int timeIdMin1,timeIdMin2;
  int timeIdMax1,timeIdMax2;

  int i;
  int i1,i2;

  /* Allocate the signals to build the two atoms in (once only) */
  if(atomR1 == NULL) {
    atomR1 = NewSignal();
    atomI1 = NewSignal();
    atomR2 = NewSignal();
    atomI2 = NewSignal();
  }

  /* 
   * Build the atoms with a unit norm in a complex signal.
   * The signals have exactly windowSize points.
   * The firstp point has a value of zero.
   */
  BuildAtom(atom1,atomR1,atomI1,borderType,YES);
  BuildAtom(atom2,atomR2,atomI2,borderType,YES);

  /* 
   * Compute the inner product :
   * a/ [timeIdMin timeIdMax] is the support
   *   (in 'absolute' sample coordinates) of the
   *   intersection of the supports of the atoms
   * b/ timeIdMin<j> is the index (in 'absolute' 
   *   sample coordinates) of the first non zero point in atom<j>
   */
  *pReal = *pImag = 0.0;
  AtomsIntersect(atom1,atom2,&timeIdMin,&timeIdMax);
  ComputeWindowSupport(atom1->windowSize,atom1->windowShape,atom1->timeId,&timeIdMin1,&timeIdMax1);
  ComputeWindowSupport(atom2->windowSize,atom2->windowShape,atom2->timeId,&timeIdMin2,&timeIdMax2);

  for(i = timeIdMin; i <= timeIdMax; i++) {
    /* 
     *Checked correct on the 22/02/2001 by R. Gribonval
     * The indexes in the atom windows : 
     * because we always have atom[0] = 0.0,
     *   i=timeIdMin1 corresponds to i1=1 
     *   i=timeIdMin2 corresponds to i2=1 
     * hence the expression of i1,i2
     * By construction we will (should?) always have
     *  1 <= i1,i2 <= atom(1|2)Size-1.
     */
    i1 = i-timeIdMin1+1;
    i2 = i-timeIdMin2+1;
      
    /* We multiply atom1[i] by Conjugate(atom2[i]) */
    *pReal += atomR1->Y[i1]*atomR2->Y[i2]+atomI1->Y[i1]*atomI2->Y[i2];
    *pImag += atomI1->Y[i1]*atomR2->Y[i2]-atomR1->Y[i1]*atomI2->Y[i2];
  }

  // If the atoms are indeed 'real' (i.e. unmodulated) 
  freqId1 = (long)atom1->freqId;
  freqId2 = (long)atom2->freqId;
  if((freqId1 == atom1->freqId) && (freqId2 == atom2->freqId)
     && (freqId1%GABOR_NYQUIST_FREQID)==0
     && (freqId2%GABOR_NYQUIST_FREQID)==0
     && atom1->chirpId == 0.0 && atom2->chirpId == 0.0)
    *pImag = 0.0;
  return;
}

/**********************************************/
/*
 * The inner product between two COMPLEX atoms.
 * The values of coeffR/coeffI/coeff2/realGG/imagGG/cosPhase/sinPhase 
 * do not influence the result.
 * - Some fast Analytic formulas are implemented in files 'atom_<window>innerprod.c'
 *  and used when possible ('gauss' and 'asym' only for the moment)
 * - Numerical evaluation is done for other shapes or when it is faster (small scales)
 * TODO : optimize the switch between analytic and numeric.
*/
/**********************************************/
static void CCAtomInnerProduct(const ATOM atom1,const ATOM atom2,char flagForceNumeric,LWFLOAT *pReal,LWFLOAT *pImag)
{
  /* Case of non-overlapping atoms */
  if(AtomsIntersect(atom1,atom2,NULL,NULL)==NO) {
    *pReal = *pImag = 0.0;
    return;
  }
  
  // Case where we want to force exact computation
  if(flagForceNumeric) {
    CCAtomNumericInnerProduct(atom1,atom2,pReal,pImag); 
    return;
  }

  /* Choosing the computation method adapted to the shape of the atoms */
  switch(atom1->windowShape) {
  case GaussWindowShape :
    // Case of two Gaussian atoms with large enough windows
    // TODO : We may have to optimize the windowSize switch
    if (atom2->windowShape == GaussWindowShape && atom1->windowSize > 16 && atom2->windowSize > 16)  {	
      CCAtomAnalyticGauss(atom1,atom2,pReal,pImag);
      return;
    }
    /* Numerical formula is better */
    CCAtomNumericInnerProduct(atom1,atom2,pReal,pImag); 
    return;    
  case FoFWindowShape :		
    // Case of two FoF atoms with large enough windows
    // TODO : We may have to optimize the windowSize switch
    /* ???? VOIR CE QUI EST PLUS RAPIDE ET JUSTE !! (PB FoF NUMERIQUE) */
    if (atom2->windowShape == FoFWindowShape && atom1->chirpId == 0.0 && atom2->chirpId == 0.0 &&
	atom1->windowSize > 32 && atom2->windowSize > 32)  {	
      CCAtomAnalyticFoF(atom1,atom2,pReal,pImag);  
      return;
    }
    /* Numerical formula */
    CCAtomNumericInnerProduct(atom1,atom2,pReal,pImag); 
    return;
  default :
    CCAtomNumericInnerProduct(atom1,atom2,pReal,pImag); 
    return;
  }
}

/* 
 * After going through this function, CheckAtom(atom) would generate an error because 
 * 'atom' does not satisfy 0<=freqId<=GABOR_NYQUIST_FREQID
 */
static void ConjugateAtom(ATOM atom) 
{
  CheckAtom(atom);
  atom->freqId = 2*GABOR_NYQUIST_FREQID-atom->freqId;
  if(atom->freqId >= 2*GABOR_NYQUIST_FREQID)
    atom->freqId -= 2*GABOR_NYQUIST_FREQID;
  atom->chirpId = -atom->chirpId;
}


/*************************************************************/
/* "Auto" inner-product of an atom with its conjugate.	     */
/*	<g_r+,g_r->					     */
/* The values of coeffR/coeffI/coeff2/realGG/imagGG/         */
/* cosPhase/sinPhase do not influence the result.            */
/*************************************************************/
void AutoAtomInnerProduct(const ATOM atom,char flagForceNumeric,LWFLOAT *pReal,LWFLOAT *pImag)
{
  static ATOM conjugate = NULL;

  CheckAtom(atom);

  conjugate = CopyAtom(atom,conjugate);
  ConjugateAtom(conjugate);
  /* Inner product of complex atom with its conjugate */
  CCAtomInnerProduct(atom,conjugate,flagForceNumeric,pReal,pImag);
  /* Fix the cases when the imaginary part is surely zero : 
   * it is the case for symmetric windows like GAUSSIAN  */
  /* TODO should simply test symmetry */
  if(atom->chirpId == 0.0 && GetFlagAsymWindowShape(atom->windowShape)==NO)
    *pImag = 0.0;
}

/*************************************************************************************/
/*
 * The inner product between two atoms.
 * (a REAL atom and a COMPLEX one)
 * The values of coeffR/coeffI/coeff2/realGG/imagGG/ do not influence the result.
 * Only the (cosPhase,sinPhase) of the first atom do.
 */
/****************************************************/
/*
 * <g_r,g_c>  =  <g'_r,g_c> / ||g'_r|| 
 *
 * where g'_r = 1/2(exp i phi <g'_r+,g_c> + exp -i phi <g'_r-,g_c>)
 * with g'_r+ and g'_r- complex atoms 
 * and  ||g'_r||^2 = 1/2 (1+1\Re{\exp{2 i \phi}<g'_r+,g'_r->})
 */
void RCAtomInnerProduct(const ATOM atomR,const ATOM atomC,char flagForceNumeric,LWFLOAT *pReal,LWFLOAT *pImag)
{
  static ATOM copyR = NULL;
  LWFLOAT cosPhase,sinPhase;
  long freqId;

  LWFLOAT x1,y1,x2,y2;
  LWFLOAT x,y;

  LWFLOAT realGG,imagGG;
  LWFLOAT cos2Phase,sin2Phase;
  
  LWFLOAT norm2;
  LWFLOAT factor;
  
  /* Checking */
  CheckAtomReal(atomR);
  CheckAtom(atomC);
  CheckTFContentCompat(atomR,atomC);

  /* Test if the supports intersect */
  if(AtomsIntersect(atomR,atomC,NULL,NULL)==NO)  {
    *pReal = *pImag = 0.0;
    return;
  }
  
  /* The inner-product before normalization 
   * <g'_r,g_c> = 1/2(exp i phi <g_r+,g_c> + exp -i phi <g_r-,g_c>) 
   */
  cosPhase = atomR->cosPhase;
  sinPhase = atomR->sinPhase;
  /* We must copy at least one of the two atoms to avoid bugs
   * in the conjugation process if atomR == atomC
   */
  copyR = CopyAtom(atomR,copyR);
  /* Inner product of the real (positive frequency) atom with the complex atom */
  CCAtomInnerProduct(copyR,atomC,flagForceNumeric,&x1,&y1);
  /* Multiply by \exp{i*phase}  */
  x = x1*cosPhase-y1*sinPhase;
  y = x1*sinPhase+y1*cosPhase;

  /* Taking the conjugate of the real atom */
  ConjugateAtom(copyR);
  /* Inner product of the real (negative frequency) atom with the complex atom */
  CCAtomInnerProduct(copyR,atomC,flagForceNumeric,&x2,&y2);
  /* Multiply by \exp{-i*phase} */
  x += x2*cosPhase+y2*sinPhase;
  y += x2*sinPhase-y2*cosPhase;

  /* Dividing by 2 */
  *pReal = .5*x;
  *pImag = .5*y;

  /* <g'_r+,g'_r-> */
  AutoAtomInnerProduct(atomR,flagForceNumeric,&realGG,&imagGG);