atom_innerprod.c

LastWave

  
  /* ||g'_r||^2 = 1/2 (1+\Re{\exp{2 i \phi}<g'_r+,g'_r->}) */
  cos2Phase	= cosPhase*cosPhase-sinPhase*sinPhase;
  sin2Phase	= 2*sinPhase*cosPhase;
  norm2 	= (1+cos2Phase*realGG-sin2Phase*imagGG)/2;
  
  /* Normalizing */
  factor	= sqrt(norm2);
  *pReal 	/= factor;
  *pImag 	/= factor;

  /* Case where the complex atom is not modulated */
  freqId = (long)atomC->freqId;
  if((freqId == atomC->freqId) && (freqId%GABOR_NYQUIST_FREQID)==0 && atomC->chirpId == 0.0) *pImag = 0.0;


}


/*************************************************************************************/
/*
 * The inner product between two REAL atoms (if not, an error is generated).
 * The values of coeffR/coeffI/coeff2/realGG/imagGG/ do not influence the result.
 * Only the (cosPhase,sinPhase) of the two atoms do.
 */
/****************************************************/
/*
 * <g_rn1,g_rn2> = 1/2(e^-i\phi2 <g_rn1,g_r2+> + e^+i\phi2 <g_rn1,g_r2->)/||g_r2|| 
 *
 * where ||g_r2||^2 = 1/2(1+\Re{\exp 2i\phi2 <g_r2+,g_r2->
 */

void RRAtomInnerProduct(const ATOM atom1,const ATOM atom2,char flagForceNumeric,LWFLOAT *pReal)
{
  static ATOM copy2 = NULL;
  LWFLOAT cosPhase,sinPhase;
  
  LWFLOAT x1,y1,x2,y2;
  
  LWFLOAT realGG,imagGG;
  LWFLOAT cos2Phase,sin2Phase;

  LWFLOAT norm2;
  LWFLOAT factor;
    
  /* Checking */
  CheckAtomReal(atom1);
  CheckAtomReal(atom2);
  CheckTFContentCompat(atom1,atom2);
  
  /* Test if the supports intersect */
  if(AtomsIntersect(atom1,atom2,NULL,NULL)==NO)  {
    *pReal = 0.0;
    return;
  }

  /* The inner-product before normalization of the second atom, i.e. betweeen
   * the first real atom and the second complex atom
   */
  cosPhase = atom2->cosPhase;
  sinPhase = atom2->sinPhase;
 /* We must copy at least one of the two atoms to avoid bugs
   * in the conjugation process if atomR == atomC
   */
  copy2 = CopyAtom(atom2,copy2);
  RCAtomInnerProduct(atom1,copy2,flagForceNumeric,&x1,&y1);
  /* Multiplying by \exp{-i\phi2}, we only care about the real part because in the end the imaginary part will be zero */
  *pReal  = x1*cosPhase+y1*sinPhase;

  /* Now the second term with the conjugate of the complex atom */
  ConjugateAtom(copy2);
  RCAtomInnerProduct(atom1,copy2,flagForceNumeric,&x2,&y2);
  /* Multiplying by \exp{+i\phi2} and adding */
  *pReal += x2*cosPhase-y2*sinPhase;
  
 /* Dividing by 2 */
  *pReal /= 2;
  
  if(*pReal == 0.0)
    return;
  
  
  /* We normalize the second real atom */
  /* ||g_r2||^2 = 1/2 (1+\Re{\exp{2 i \phi2}<g_r2+,g_r2->}) */
  
  /* we use atom2 and not copy2 because we need
   * <g_r2+,g_r2-> and not <g_r2-,g_r2+>
   * (copy2 is now the conjugate of atom2) 
   */
  AutoAtomInnerProduct(atom2,flagForceNumeric,&realGG,&imagGG);
  
  cos2Phase	= cosPhase*cosPhase-sinPhase*sinPhase;
  sin2Phase	= 2*sinPhase*cosPhase;
  norm2 	= (1+cos2Phase*realGG-sin2Phase*imagGG)/2;
  
  /* Normalizing */
  factor	= sqrt(norm2);
  *pReal 	/= factor;
}

/*************************************************************************************/
/*
 * The inner product between a SIGNAL and a COMPLEX atom
 * (if not, an error is generated).
 * The values of coeffR/coeffI/coeff2/realGG/imagGG/cosPhase/sinPhase
 * of the atom do not influence the result.
 */
// If the frequency is either 0 or Nyquist the atom takes indeed real
// values and we garantee that the imaginary part of the inner product is zero.
/****************************************************/
void SCAtomInnerProduct(const SIGNAL signal,const ATOM atom,char borderType,LWFLOAT *pReal,LWFLOAT *pImag)
{
  LWFLOAT real,imag;

  static SIGNAL atomR = NULL;
  static SIGNAL atomI = NULL;
  static SIGNAL signalExtract = NULL;

  int timeIdMin,timeIdMax,i;

  // Checking 
  if(signal == NULL) Errorf("SCAtomInnerProduct : NULL signal");
  CheckAtom(atom);
  if(!BorderTypeIsOK(borderType)) Errorf("SCAtomInnerProduct : bad borderType '%d'",borderType);

  // Initialization
  real = imag = 0.0;
  /* Are the supports of the signal and the atom overlapping ? */
  ComputeWindowSupport(atom->windowSize,atom->windowShape,atom->timeId,&timeIdMin,&timeIdMax);
  switch(borderType) {
  case BorderPad0 :      
    if ((timeIdMax < 0) || (timeIdMin > signal->size-1))
      return;
    break;
  default : 
    Errorf("SCAtomInnerProduct : borderType '%s' not treated yet",BorderType2Name(borderType));
    break;
  }

  /* Allocation (just once) */
  if(atomR == NULL) {
    atomR = NewSignal();
    atomI = NewSignal();
    signalExtract = NewSignal();
  }

  /* 
   * Build the atom with a unit norm in a complex signal.
   * The signal has exactly windowSize points.
   * The firstp point has a value of zero.
   */
  BuildAtom(atom,atomR,atomI,borderType,YES);

  /* 
   * Extract the corresponding piece of 'signal'
   * The first point must coincide with the first point
   * of the window of the atom, which is timeIdMin-1.
   */
  ExtractSig(signal,signalExtract,borderType,timeIdMin-1,atomR->size);

  /* Compute the real part of the inner product */
  for(i=0; i< atomR->size; i++) 
    real += signalExtract->Y[i]*atomR->Y[i];
  /* Compute the imaginary part if necessary (i.e. if nonzero) */
  if(atom->freqId != 0 && atom->freqId != GABOR_NYQUIST_FREQID) 
    for(i=0; i< atomR->size; i++) 
      imag -= signalExtract->Y[i]*atomI->Y[i];

  *pReal=real;
  *pImag=imag;
}


/*****************************************************************************/
/*                            TESTS  !!!!                                    */
/*****************************************************************************/

void C_Inner(char **argv)
{
    ATOM atom1=NULL,atom2=NULL;
    SIGNAL signal=NULL;
    char flagForceNumeric = NO;
    char *action;
    char *name;
    char opt;

    LISTV lv;
    /* Result */
    int borderType = STFT_DEFAULT_BORDERTYPE;
    LWFLOAT re,im;

    /* Parsing */
    argv = ParseArgv(argv,tWORD,&action,-1);

    /*
     * Complex-signal
     */    
    if(!strcmp(action,"sig")) {
      argv = ParseArgv(argv,tSIGNAL,&signal,tATOM,&atom1,-1);
      while( (opt = ParseOption(&argv))) {
	switch(opt) {
	case 'b' :
	  argv = ParseArgv(argv,tSTR,&name,-1);
	  borderType = Name2BorderType(name);
	  break;
	default :
	  ErrorOption(opt);
	}
      }
      NoMoreArgs(argv);
      SCAtomInnerProduct(signal,atom1,borderType,&re,&im);
      lv = TNewListv();
      AppendFloat2Listv(lv,re);
      AppendFloat2Listv(lv,im);
      SetResultValue(lv);
      return;
    }

    /*
     * "Auto" inner-product with the conjugate 
     */
    if(!strcmp(action,"auto")) {
      argv = ParseArgv(argv,tATOM,&atom1,0);
      AutoAtomInnerProduct(atom1,NO,&re,&im);
      lv = TNewListv();
      AppendFloat2Listv(lv,re);
      AppendFloat2Listv(lv,im);
      SetResultValue(lv);
      return;
    }
    
    /*
     * Inner product between two atoms 
     */
    argv = ParseArgv(argv,tATOM,&atom1,tATOM,&atom2,-1);
    while( (opt = ParseOption(&argv))) {
      switch(opt) {
      case 'n' :
	flagForceNumeric = YES;
	NoMoreArgs(argv);
	break;
      default :
	ErrorOption(opt);
      }
    }
    /*
     * Complex-Complex
     */
    if(!strcmp(action,"cc"))  {
      CheckTFContentCompat(atom1,atom2);
      CCAtomInnerProduct(atom1,atom2,flagForceNumeric,&re,&im);
      lv = TNewListv();
      AppendFloat2Listv(lv,re);
      AppendFloat2Listv(lv,im);
      SetResultValue(lv);
      return;
    }
    /* 
     * Real-complex 
     */
    if(!strcmp(action,"rc")) {
      RCAtomInnerProduct(atom1,atom2,flagForceNumeric,&re,&im);
      lv = TNewListv();
      AppendFloat2Listv(lv,re);
      AppendFloat2Listv(lv,im);
      SetResultValue(lv);
      return;
    }
    /*
     * Real-real normalized
     */
    if(!strcmp(action,"rr")) {
      RRAtomInnerProduct(atom1,atom2,flagForceNumeric,&re);
      SetResultFloat(re);
      return;
    }
    Printf("Unknown action '%s'",action);
    ErrorUsage();
}

/* EOF */