pf_lib.c

来自「LastWave」· C语言 代码 · 共 1,800 行 · 第 1/4 页

    flagSwap = PFYes;
  else
    flagSwap = PFNo;
    
  /* If pf->cellArray != NULL it means that we are writing
     on an already used PartitionFunction
     So we delete the cellArray */
  if( pf->cellArray != NULL )
    {
      PFCellArrayDelete(pf->cellArray,pf->qNumber);
      pf->cellArray = NULL;
    }
  
  
  if( PFReadBinHeader(fp,pf,flagSwap) != PFYes)
    return PFReadErr;
  
  /* We verify that the values are meaningful */
  if( _PFIsValid(pf)  != PFYes )
    return PFNotValid;

  /* We allocate space for the qArray */
  qArray = (double *) calloc(pf->qNumber,sizeof(double));
  if( qArray == NULL ) return PFErrAlloc;

  
  /* We initialize the cell array */
  pf->cellArray = PFCellArrayNew(pf->octaveNumber*pf->voiceNumber,
				 pf->qNumber,qArray);
  if( pf->cellArray == NULL) 
    {
      free(qArray);
      return PFErrAlloc;
    }
  
  /* We read the content of the cells */
  for(nq=0;nq < pf->qNumber;nq++)
    if( PFCellReadBin(fp,pf->cellArray[nq],flagSwap) != PFYes )
      {
	free(qArray);
	return PFReadErr;
      }

  free(qArray);
  return PFYes;
}

int PFRead(FILE *fp,PartitionFunction pf)
{
  char tempLine[PFMESSMAXSIZE+2];

  if(pf == NULL) return PFNo;
  
  if(fp == NULL) return PFNo;
  
  /* Are we dealing with a partition function file ? */
  if( fgets(tempLine,PFMESSMAXSIZE+2,fp) == NULL)
    return PFReadErr;
  /* The last character in templine is \n we convert it to \0
     before comparing strings */
  tempLine[strlen(tempLine)-1] = '\0';
  if( strcmp(tempLine,PFMessIdent) != 0 )
    return PFErrFormat;

  /* Is it an ASCII or a BINARY file ? */
  if( fgets(tempLine,PFMESSMAXSIZE+2,fp) == NULL)
    return PFReadErr;
  /* The last character in templine is \n we convert it to \0
     before comparing strings */
  tempLine[strlen(tempLine)-1] = '\0';
  if( strcmp(tempLine,PFMessAsciiFormat) == 0 )
    return PFReadAscii(fp,pf);
  else if( strcmp(tempLine,PFMessBinFormat) == 0 )
    return PFReadBin(fp,pf);
  else return PFErrFormat;
}


/*************************************************
 *
 * computation of the partition functions at one scale
 *
 *************************************************/

/* The ordering of data in PFComputeOneScaleF() is done using this function.
   One can change it using PFChangeQsort() especially if no sorting is wanted.
   The new sorting function must work the same as qsort() of <stdlib.h>
   */
static void PFDefaultQsortData(double *array,int size)
{
  qsort((void *)array,size,sizeof(double),
	(int (*)(const void *,const void *))PFCompDouble);
}
void PFMinimalQsortData(double *array,int size)
{
  /* Minimal sorting consist in :
     - putting all 0 elt first
     - the minimal value next
     - the maximal value as last element
     We assume all elt are >= 0
  */
  int i;
  int imin = 0;
  int imax = size-1;
  double temp;
#define PFSWAP(a,b) temp = (a); (a) = (b); (b) = temp

  /**** initialisation */
  /* Le's be sure that max is not 0 */
  i = imax;
  while(i >= 0 && array[i] == 0)
    i--;
  if(i < 0)
    return;
  PFSWAP(array[i],array[imax]);
  /* Le's be sure that we can find a none 0 element at imin */
  while(imin < size && array[imin] == 0)
    imin++;
  if(imin >= size-1)
    return;
  /* min should be smaller than max */
  if(array[imin] > array[imax]) {
    PFSWAP(array[imin],array[imax]);
  }
  
  for(i=imin+1; i<imax; i++) {
    if(array[i] == 0) {
      array[i] = array[imin+1];
      array[imin+1] = array[imin];
      array[imin] = 0;
      imin++;
    } else if (array[i] < array[imin]) {
      PFSWAP(array[i],array[imin]);
    } else if (array[i] > array[imax]) {
      PFSWAP(array[i],array[imax]);
    }
  }
}
#undef PFSWAP
static void (*PFQsortData)(double *array,int  size) = PFDefaultQsortData;

void PFChangeQsort( void (*NewQsort)(double *,int ))
{
  if(NewQsort == NULL)
    PFQsortData = PFDefaultQsortData;
  else
    PFQsortData = NewQsort;
}

/* LWFLOAT version: the wavelet coefficients are given 
   in an array of LWFLOAT */
/* scales go from 0 to octaveNumber*voiceNumber-1 */
int PFComputeOneScaleFLOAT(PartitionFunction pf,int scale,const LWFLOAT *t,int tSize)
{
  int nq,i,imin;
  double tm;
  double q;
  double *tempAbsT,*tempTq,*tempLogT;
  PartitionFunctionCell tempPfc;

  /* Some verifications */
  if(pf == NULL) 
    return PFNo;
  if( (PFIsValid(pf) != PFYes) || (pf->qNumber <= 0) || (scale < 0) 
      || (scale >= pf->octaveNumber*pf->voiceNumber) || (tSize <= 0) )
    return PFNotValid;
  
  /* One is supposed to compute scales in increasing order */
  if( (scale - pf->indexMax) != 1 )
    return PFNotValid;

  /* Allocation of space for tempAbsT and tempMonome
     using only one call to malloc */
  tempAbsT = (double *) malloc(3*tSize*sizeof(double));
  if(tempAbsT == NULL) return PFNo;
  tempTq = tempAbsT + tSize;
  tempLogT = tempAbsT + 2*tSize;

  /* Computation of |t|, sorting of |t| and determination of imin 
     such that for i >= imin |t|(i) != 0 
     */
  for(i=0;i<tSize;i++)
    tempAbsT[i] = fabs((double) t[i]);
  PFQsortData(tempAbsT,tSize);
  imin = 0;
  while( imin < tSize && tempAbsT[imin] == 0 )
    imin++;

  /*  Printf("imin= %d, Tmin= %g, Tmax %g\n",imin,tempAbsT[imin],
      tempAbsT[tSize-1]);
      */

  /* If imin = tSize then all Ts are equal to zero 
     => nothing to compute
     (Initialisation to 0 of PFs is done at alloc with calloc() ) */

  if (imin < tSize) {
    for(nq=0; nq < pf->qNumber; nq++)
      {
	tempPfc = pf->cellArray[nq];
	q = tempPfc->q;
	
	/* Do we want T/tm to be >= 1 or <= 1 */
	if(q >= 0.)
	  tm = tempAbsT[imin];
	else
	  tm = tempAbsT[tSize-1];
	
	/* We compute Tq and Log(T/tm) */
	for(i=imin; i < tSize; i++)
	  {
	    tempTq[i] = pow(tempAbsT[i],q);
	    tempLogT[i] = log(tempAbsT[i]/tm);
	  }
	
	/* We compute sTq and sTqLogT */
	if(q >= 0.)
	  {
	    for(i=imin; i < tSize; i++)
	      {
		/* The arrays of a PartitionFunctionCell are obtained with
		   calloc, so there is no need to initialize them
		*/
		tempPfc->sTq[scale] += tempTq[i];
		tempPfc->sTqLogT[scale] += tempTq[i]*tempLogT[i];
	      }
	  }
	else
	  {
	    for(i=tSize-1; i >= imin; i--)
	      {
		/* The arrays of a PartitionFunctionCell are obtained with
		   calloc, so there is no need to initialize them
		*/
		tempPfc->sTq[scale] += tempTq[i];
		tempPfc->sTqLogT[scale] += tempTq[i]*tempLogT[i];
	      }
	  }
	/* sTqLogT = sTqLog(T/tm)+log(tm)*sTq
	 */
	tempPfc->sTqLogT[scale] = 
	  tempPfc->sTqLogT[scale] + log(tm)*tempPfc->sTq[scale];
	
	/* We compute LogSTq, sTqLogT_sTq, log2STq, sTqLogT_sTq2 
	   and logSTqSTqLogT_sTq
	*/
	if(tempPfc->sTq[scale] != 0.0)
	  {
	    tempPfc->logSTq[scale] = 
	      log(tempPfc->sTq[scale]/((double) (tSize-imin)));
	    tempPfc->log2STq[scale] = 
	      tempPfc->logSTq[scale]*tempPfc->logSTq[scale];
	    
	    tempPfc->sTqLogT_sTq[scale] = 
	      tempPfc->sTqLogT[scale]/tempPfc->sTq[scale];
	    tempPfc->sTqLogT_sTq2[scale] = 
	      tempPfc->sTqLogT_sTq[scale]*tempPfc->sTqLogT_sTq[scale];
	    tempPfc->logSTqSTqLogT_sTq[scale] = 
	      tempPfc->logSTq[scale]*tempPfc->sTqLogT_sTq[scale];
	  }
	else 
	  {
	    tempPfc->logSTq[scale] = 0.0;
	    tempPfc->sTqLogT_sTq[scale] = 0.0;
	    tempPfc->log2STq[scale] = 0.0;
	    tempPfc->sTqLogT_sTq2[scale] = 0.0;
	    tempPfc->logSTqSTqLogT_sTq[scale] = 0.0;
	  }
	
      }
  }

  pf->indexMax = scale;
  free(tempAbsT);
  return PFYes;
}


/*************************************************
 *
 * Addition of two PartitionFunctions
 *
 *************************************************/

/* Comparaison of qList */
enum PFCompQlistCode { PFQListEqual, PFQListDifferent, PFQListMixed };
static int PFCompQList(const PartitionFunction pf1,const PartitionFunction pf2)
{
  int nq1,nq2;
  int flagEqual = PFNo;

  assert( PFIsValid(pf1) == PFYes );
  assert( PFIsValid(pf2) == PFYes );
  
  /* Are they strickly equal ? */
  if( pf1->qNumber == pf2->qNumber)
    {
      flagEqual = PFYes;
      for(nq1=0;nq1 < pf1->qNumber;nq1++)
	if(fabs(pf1->cellArray[nq1]->q - pf2->cellArray[nq1]->q) > PFQDIFFMAX)
	   {
	     flagEqual = PFNo;
	     break;
	   }
    }
  if(flagEqual == PFYes)
    return PFQListEqual;

  /* Are they strickly different ? */
  for(nq1=0;nq1 < pf1->qNumber;nq1++)
      for(nq2=0;nq2 < pf2->qNumber;nq2++)
	{
	  if(fabs(pf1->cellArray[nq1]->q - pf2->cellArray[nq2]->q) <= PFQDIFFMAX)
	    return PFQListMixed;
	}

  return PFQListDifferent;
}

/* Addition of two PartitionFunctionCell, the result is in pfc1.
   q's should be the same but pfc1->size may be bigger.
   ( it is assumed that the difference only comes from a different
   octaveNumber ) 
   */
static int PFCellAdd(PartitionFunctionCell pfc1,
		     const PartitionFunctionCell pfc2)
{
  int i;

  if( (pfc1->size < pfc2->size) || fabs(pfc1->q - pfc2->q) > PFQDIFFMAX ) 
    return PFNotCompatible;

  for(i=0;i < pfc2->size;i++)
    {
      pfc1->sTq[i] += pfc2->sTq[i];
      pfc1->sTqLogT[i] += pfc2->sTqLogT[i];
      pfc1->logSTq[i] += pfc2->logSTq[i];
      pfc1->sTqLogT_sTq[i] += pfc2->sTqLogT_sTq[i];
      pfc1->log2STq[i] += pfc2->log2STq[i];
      pfc1->sTqLogT_sTq2[i] += pfc2->sTqLogT_sTq2[i];
      pfc1->logSTqSTqLogT_sTq[i] += pfc2->logSTqSTqLogT_sTq[i];
      
    }    
  
  return PFYes;
}
/* The result will be in pf1 */
int PFStandardAddition(PartitionFunction pf1,const PartitionFunction pf2)
{
  int nq;
  int compQList;
  PartitionFunctionCell *cellArray;

  if( (PFIsValid(pf1) != PFYes) || (PFIsValid(pf2) != PFYes) )
    return PFNotValid;

  /* method, aMin, octaveNumber, voiceNumber, signalSize, dimension
     have to be equal for the addition to be meaningful */
  if( (strcmp(pf1->method,pf2->method) != 0) 
      || (fabs(pf1->aMin - pf2->aMin) > PFQDIFFMAX)
      || (pf1->octaveNumber != pf2->octaveNumber)
      || (pf1->voiceNumber != pf2->voiceNumber) 
      || (pf1->signalSize != pf2->signalSize) 
      || (pf1->dimension != pf2->dimension) )
    return PFNotCompatible;
  
  compQList = PFCompQList(pf1,pf2);

  /* If qList are mixed we can't add */
  if( compQList == PFQListMixed )
    return PFNotCompatible;
  /* If the qList are strickly different and all the other fields are 
     the same (except indexMax and qNumber) then the result is simply 
     the concatenation of the two cellArrays: we are adding new values of q
     */
  else if( compQList == PFQListDifferent )
    {
      if(pf1->signalNumber != pf2->signalNumber)
	return PFNotCompatible;
      
      cellArray = PFCellArrayConcat(pf1->qNumber,pf1->cellArray,
				    pf2->qNumber,pf2->cellArray);
      if(cellArray == NULL) return PFErrAlloc;

      PFCellArrayDelete(pf1->cellArray,pf1->qNumber);
      pf1->cellArray = cellArray;
      pf1->qNumber = pf1->qNumber + pf2->qNumber;
     
      if(pf1->indexMax > pf2->indexMax)
	pf1->indexMax = pf2->indexMax;
      
      return PFYes;
    }
  /* If the qList are strickly equal and all the other fields are 
     the same (except indexMax and signalNumber) then the result is
     the addition of the partition function for each q: we are adding 
     more statistic 
     */
  else if( compQList == PFQListEqual )
    {
      for(nq=0;nq < pf1->qNumber;nq++)
	if( PFCellAdd(pf1->cellArray[nq],pf2->cellArray[nq]) != PFYes)
	  return PFNotCompatible;
      pf1->signalNumber = pf1->signalNumber + pf2->signalNumber;

      if(pf1->indexMax > pf2->indexMax)
	pf1->indexMax = pf2->indexMax;

      return PFYes;
    }

  return PFNo;
}

int PFNonStandardAddition(PartitionFunction pf1,const PartitionFunction pf2)
{
  int compQList;
  int nq;
  
  if( (PFIsValid(pf1) != PFYes) || (PFIsValid(pf2) != PFYes) )
    return PFNotValid;

  compQList = PFCompQList(pf1,pf2);
  
  /* method, aMin, voiceNumber, dimension, qList must be equal for
     the addition to be meaningful.
     octaveNumber2 must be  smaller than octaveNumber1 */
  if( (strcmp(pf1->method,pf2->method) != 0) 
      || (fabs(pf1->aMin - pf2->aMin) > PFQDIFFMAX)
      || (pf1->octaveNumber < pf2->octaveNumber)
      || (pf1->voiceNumber != pf2->voiceNumber)
      || (pf1->dimension != pf2->dimension)
      || (compQList != PFQListEqual)
      )
    return PFNotCompatible;

  for(nq=0;nq < pf1->qNumber;nq++)
    if( PFCellAdd(pf1->cellArray[nq],pf2->cellArray[nq]) != PFYes)
      return PFNotCompatible;

  pf1->signalSize = pf1->signalSize + pf2->signalSize;
  pf1->signalNumber = pf1->signalNumber + pf2->signalNumber;
  

  if(pf1->indexMax > pf2->indexMax)
    pf1->indexMax = pf2->indexMax;
  
  return PFYes;
}