/*..........................................................................*//*                                                                          *//*      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  : remi.gribonval@inria.fr                                    *//*               lastwave@cmapx.polytechnique.fr                            *//*                                                                          *//*..........................................................................*//*                                                                          *//*      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,