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



/****************************************************/
/*
 * 	The HIGHRES Short Time Fourier Transform :
 *
 *		COMPUTATION and UPDATE
 */
/****************************************************/
#include "lastwave.h"

#include "atom.h"

/*#define DEBUG_HIGHRES*/

/*************************************************************/
/*
 *	COMPUTATION OF THE HIGH-RESOLUTION CORRELATION
 */
/*************************************************************/

extern void RRAtomInnerProduct(const ATOM atom1,const ATOM atom2,char flagForceNumeric,LWFLOAT *pReal);

/******************************/
/* 
 *	The function to compute 
 * 	the high-res correlation
 */
/******************************/
/*
 * We assume that the atom ENERGY and PHASE are already properly set
 * [firstp,lastp] is the range of timeId not affected by border effects
 * in the HIGHRES stft at the atom's windowSize
 */
static void SetAtomCoeff2HighRes(ATOM atom,STFT stftHighRes,STFT subStftReal,STFT subStftPhase)
{
  int timeId;
  
  static ATOM subAtom = NULL;
  int timeIdMinAtom,timeIdMaxAtom;
  int shiftMin,shiftMax;
  int timeIdMinSubAtoms,timeIdMaxSubAtoms;
  int tmp;

  int subTimeIdMin,subTimeIdMax;

  int subTRate;
  

  int subTimeId;
  int subFreqId;
  
  LWFLOAT *subCoeff2s,*subPhases;
  LWFLOAT subCoeff2,subPhase;
  LWFLOAT crossInner;
  
  /* Checking */
  CheckAtomReal(atom);
  if(!BorderTypeIsOK(stftHighRes->borderType))
    Errorf("SetAtomCoeff2HighRes : bad borderType %d",stftHighRes->borderType);
  CheckStftReal(subStftReal);
  CheckStftPhase(subStftPhase);
  if(subStftReal->flagUpToDate == NO || subStftPhase->flagUpToDate == NO)
    Errorf("SetAtomCoeff2HighRes: subStftReal or subStftPhase is out of date");
  // The atom MUST be on the time-frequency grid of subStft 
  // because we cannot recompute the sub-atom's energy
  timeId = (int) atom->timeId;
  if(atom->timeId != timeId || timeId % subStftReal->tRate != 0)
    Errorf("SetAtomCoeff2HighRes: timeId %g not on grid %d",atom->timeId,subStftReal->tRate);
  if(atom->freqId != (int)atom->freqId || ((int)atom->freqId) % subStftReal->fRate != 0)
    Errorf("SetAtomCoeff2HighRes: freqId %g not on grid %d",atom->freqId,subStftReal->fRate);
  
  /* Allocating (once) */
  if(subAtom == NULL) {
    subAtom = NewAtom();
  }
  else {
    ClearAtom(subAtom);
  }
  CopyFieldsTFContent(subStftReal,subAtom);
  subAtom->windowShape  = subStftReal->windowShape;
  subAtom->windowSize   = subStftReal->windowSize;
	subAtom->freqId    = (int) atom->freqId;
  subFreqId  =  (int) atom->freqId;
  /* Case when the energy is ZERO */
  if (atom->coeff2 == 0.0) return;
  
  // The support of the main atom is [timeIdMinAtom,timeIdMaxAtom]
  ComputeWindowSupport(atom->windowSize,atom->windowShape,atom->timeId,
		     &timeIdMinAtom,&timeIdMaxAtom);
  // A subAtom at 'subTimeId' has support 
  //[subTimeId+shiftMin,subTimeId+shiftMax]
  ComputeWindowSupport(subAtom->windowSize,subAtom->windowShape,0,
		     &shiftMin,&shiftMax);
  // This support is included in [timeIdMinAtom,timeIdMaxAtom] iff
  // subTimeId+shiftMin >= timeIdMinAtom
  //    and
  // subTimeId+shiftMax <= timeIdMaxAtom
  subTimeIdMin = timeIdMinAtom-shiftMin;
  subTimeIdMax = timeIdMaxAtom-shiftMax;

  // SubAtoms can only be on the grid at rate subTRate
  subTRate     = subAtom->windowSize;  //    subTRate = subStftReal->tRate;

  // Treat BORDER EFFECT : 
  // which subAtoms shall we actually look for
  switch(stftHighRes->borderType) {
    case BorderPad0 :
      // The subAtoms ??? CHANGE TODO ALL THIS BORDER EFFECT PART
      QuantizeRangeLarge(subTimeIdMin,subTimeIdMax,subTRate,
			 &subTimeIdMin,&subTimeIdMax);
      subTimeIdMin = MAX(subTimeIdMin,0);
      subTimeIdMax = MIN(subTimeIdMax,(subStftReal->nFrames*subTRate)-subTRate);

      // The union of the supports of the subAtoms is
      // [timeIdMinSubAtoms,timeIdMaxSubAtoms]
      ComputeWindowSupport(subAtom->windowSize,subAtom->windowShape,subTimeIdMin,
			 &timeIdMinSubAtoms,&tmp);
      ComputeWindowSupport(subAtom->windowSize,subAtom->windowShape,subTimeIdMax,
			 &tmp,&timeIdMaxSubAtoms);
      
      // The support [timeIdMinAtom,timeIdMaxAtom] 
      // of the atom must be ENTIRELY in [0 signalSize-1]
      if ((timeIdMinAtom < 0) || (timeIdMaxAtom > stftHighRes->signalSize-1)) {
	atom->coeff2 = 0.0;
	return;
      }
      // The support of all sub-atoms should satisfy the same
      if((timeIdMinSubAtoms < 0) ||
	 (timeIdMaxSubAtoms > subStftReal->signalSize-1)) {
	atom->coeff2 = 0.0;
	return;
      }
      break;
  default :
    Errorf("SetAtomCoeff2Highres : border type %s not treated yet",BorderType2Name(stftHighRes->borderType));
  }

  // Loop on the subAtoms timeIds : 
  // -initially we set atom->coeff2 to its maximum possible value
  // -we decrease it if necessary
  for(subTimeId = subTimeIdMin; 
      subTimeId <= subTimeIdMax; 
      subTimeId += subTRate) {
    /* The coefficient of the subAtom */
    GetStftData(subStftReal,subTimeId,&subCoeff2s,NULL);
    GetStftData(subStftPhase,subTimeId,&subPhases,NULL);
    subCoeff2 = subCoeff2s[subFreqId/subStftReal->fRate];
    subPhase  = subPhases[subFreqId/subStftPhase->fRate];
    
    // Properly setting the remaining fields of the sub-atom
    // to compute its
    // inner-product with the "big" one
    subAtom->timeId = subTimeId;
    if(!ComputeWindowGG(subAtom->windowShape,subAtom->windowSize,subAtom->freqId,&(subAtom->realGG),&(subAtom->imagGG)))
      Errorf("SetAtomCoeff2HighRes : error with ComputeWindowGG");

    subAtom->coeff2 = subCoeff2;
    // TODO : Is it costly ?? Do we need it ? Should we try something cheaper ?? 
    subAtom->cosPhase= cos(2*M_PI*subPhase);
    subAtom->sinPhase= sin(2*M_PI*subPhase);
    RRAtomInnerProduct(atom,subAtom,NO,&crossInner);

    // When the subAtom is orthogonal to the large one
    // it is NOT taken into account.
    if(crossInner == 0.0) {
      continue;
    }

    // If the sub-atom coefficient is zero !
    if(subCoeff2 == 0.0) {
      atom->coeff2 = 0.0;
      return;
    }

    // Case of a change of sign/phase
    if(crossInner < 0.0) {
      atom->coeff2 = 0.0;
      return;
    }

    // Taking the minimum encountered value
    if(atom->coeff2 > subCoeff2/(crossInner*crossInner))
      atom->coeff2 = subCoeff2/(crossInner*crossInner);
  }
}

/*************************************************/
/*
 *           Update of the REAL HIGHRES stft data,
 *           given a REAL ENERGY stft. 
 *           We loop on time, within a SPECIFIED TIME RANGE
 */
/*************************************************/

static void CheckStftHighResEnergy(STFT stftHighRes,