tns3.c

MPEG2/MPEG4编解码参考程序(实现了MPEG4的部分功能)

    }
  }  
}

/*****************************************************/
/* TnsFilter:                                        */
/*   Filter the given spec with specified length     */
/*   using the coefficients specified in filter.     */
/*   Not that the order and direction are specified  */
/*   withing the TNS_FILTER_DATA structure.          */
/*****************************************************/
void TnsFilter(int length,double* spec,TNS_FILTER_DATA* filter) {
	int i,j,k=0;
	int order=filter->order;
	double* a=filter->aCoeffs;

	/* Determine loop parameters for given direction */
	if (filter->direction) {

		/* Startup, initial state is zero */
		for (i=length-2;i>(length-1-order);i--) {
			k++;
			for (j=1;j<=k;j++) {
				spec[i]-=spec[i+j]*a[j];
			}
		}
		
		/* Now filter completely inplace */
		for	(i=length-1-order;i>=0;i--) {
			for (j=1;j<=order;j++) {
				spec[i]-=spec[i+j]*a[j];
			}
		}


	} else {

		/* Startup, initial state is zero */
		for (i=1;i<order;i++) {
			for (j=1;j<=i;j++) {
				spec[i]-=spec[i-j]*a[j];
			}
		}
		
		/* Now filter completely inplace */
		for	(i=order;i<length;i++) {
			for (j=1;j<=order;j++) {
				spec[i]-=spec[i-j]*a[j];
			}
		}
	}
}


/********************************************************/
/* TnsInvFilter:                                        */
/*   Inverse filter the given spec with specified       */
/*   length using the coefficients specified in filter. */
/*   Not that the order and direction are specified     */
/*   withing the TNS_FILTER_DATA structure.             */
/********************************************************/
void TnsInvFilter(int length,double* spec,TNS_FILTER_DATA* filter) {
	int i,j,k=0;
	int order=filter->order;
	double* a=filter->aCoeffs;
	double* temp;

    temp = (double *) malloc (length * sizeof (double));
    if (!temp) {
      CommonExit( 1, "TnsInvFilter: Could not allocate memory for TNS array\nExiting program...\n");
    }

	/* Determine loop parameters for given direction */
	if (filter->direction) {

		/* Startup, initial state is zero */
		temp[length-1]=spec[length-1];
		for (i=length-2;i>(length-1-order);i--) {
			temp[i]=spec[i];
			k++;
			for (j=1;j<=k;j++) {
				spec[i]+=temp[i+j]*a[j];
			}
		}
		
		/* Now filter the rest */
		for	(i=length-1-order;i>=0;i--) {
			temp[i]=spec[i];
			for (j=1;j<=order;j++) {
				spec[i]+=temp[i+j]*a[j];
			}
		}


	} else {

		/* Startup, initial state is zero */
		temp[0]=spec[0];
		for (i=1;i<order;i++) {
			temp[i]=spec[i];
			for (j=1;j<=i;j++) {
				spec[i]+=temp[i-j]*a[j];
			}
		}
		
		/* Now filter the rest */
		for	(i=order;i<length;i++) {
			temp[i]=spec[i];
			for (j=1;j<=order;j++) {
				spec[i]+=temp[i-j]*a[j];
			}
		}
	}
	free(temp);
}





/*****************************************************/
/* TruncateCoeffs:                                   */
/*   Truncate the given reflection coeffs by zeroing */
/*   coefficients in the tail with absolute value    */
/*   less than the specified threshold.  Return the  */
/*   truncated filter order.                         */
/*****************************************************/
int TruncateCoeffs(int fOrder,double threshold,double* kArray) {
  int i;

  for (i=fOrder;i>=0;i--) {
    kArray[i] = (fabs(kArray[i])>threshold) ? kArray[i] : 0.0;
    if (kArray[i]!=0.0) {
      return i;
    }
  }
  return ( 0 );
}

/*****************************************************/
/* QuantizeReflectionCoeffs:                         */
/*   Quantize the given array of reflection coeffs   */
/*   to the specified resolution in bits.            */
/*****************************************************/
void QuantizeReflectionCoeffs(int fOrder,
			      int coeffRes,
			      double* kArray,
			      int* indexArray) {

	double iqfac,iqfac_m;
	int i;

	iqfac = ((1<<(coeffRes-1))-0.5)/(PI/2);
	iqfac_m = ((1<<(coeffRes-1))+0.5)/(PI/2);

	/* Quantize and inverse quantize */
	for (i=1;i<=fOrder;i++) {
		indexArray[i] = (int)(0.5+(asin(kArray[i])*((kArray[i]>=0)?iqfac:iqfac_m)));
		kArray[i] = sin((double)indexArray[i]/((indexArray[i]>=0)?iqfac:iqfac_m));
	}
}

/*****************************************************/
/* Autocorrelation,                                  */
/*   Compute the autocorrelation function            */
/*   estimate for the given data.                    */
/*****************************************************/
void Autocorrelation(int maxOrder,        /* Maximum autocorr order */
		     int dataSize,		  /* Size of the data array */
		     double* data,		  /* Data array */
		     double* rArray) {	  /* Autocorrelation array */

  int order,index;

  for (order=0;order<=maxOrder;order++) {
    rArray[order]=0.0;
    for (index=0;index<dataSize;index++) {
      rArray[order]+=data[index]*data[index+order];
    }
    dataSize--;
  }
}



/*****************************************************/
/* LevinsonDurbin:                                   */
/*   Compute the reflection coefficients for the     */
/*   given data using LevinsonDurbin recursion.      */
/*   Return the prediction gain.                     */
/*****************************************************/
double LevinsonDurbin(int fOrder,          /* Filter order */
					  int dataSize,		   /* Size of the data array */
					  double* data,		   /* Data array */
					  double* kArray) {	   /* Reflection coeff array */

	int order,i;
	double signal;
	double error, kTemp;				/* Prediction error */
	double aArray1[TNS_MAX_ORDER+1];	/* Predictor coeff array */
	double aArray2[TNS_MAX_ORDER+1];	/* Predictor coeff array 2 */
	double rArray[TNS_MAX_ORDER+1];		/* Autocorrelation coeffs */
	double* aPtr = aArray1;				/* Ptr to aArray1 */
	double* aLastPtr = aArray2;			/* Ptr to aArray2 */
	double* aTemp;

	/* Compute autocorrelation coefficients */
	Autocorrelation(fOrder,dataSize,data,rArray);
	signal=rArray[0];	/* signal energy */

	/* Set up pointers to current and last iteration */ 
	/* predictor coefficients.						 */
	aPtr = aArray1;
	aLastPtr = aArray2;
	/* If there is no signal energy, return */
	if (!signal) {
	  kArray[0]=1.0;
	  for (order=1;order<=fOrder;order++) {
	    kArray[0]=0.0;
	  }
	  return 0;
	
	} else {

	  /* Set up first iteration */
	  kArray[0]=1.0;
	  aPtr[0]=1.0;		/* Ptr to predictor coeffs, current iteration*/
	  aLastPtr[0]=1.0;	/* Ptr to predictor coeffs, last iteration */
	  error=rArray[0];
	  
	  /* Now perform recursion */
	  for (order=1;order<=fOrder;order++) {
	    kTemp = aLastPtr[0]*rArray[order-0];
	    for (i=1;i<order;i++) {
	      kTemp += aLastPtr[i]*rArray[order-i];
	    }
	    kTemp = -kTemp/error;
	    kArray[order]=kTemp;
	    aPtr[order]=kTemp;
	    for (i=1;i<order;i++) {
	      aPtr[i] = aLastPtr[i] + kTemp*aLastPtr[order-i];
	    }
	    error = error * (1 - kTemp*kTemp);
	    
	    /* Now make current iteration the last one */
	    aTemp=aLastPtr;
	    aLastPtr=aPtr;		/* Current becomes last */
	    aPtr=aTemp;			/* Last becomes current */
	  }
	  return signal/error;	/* return the gain */
	}
}


/*****************************************************/
/* StepUp:                                           */
/*   Convert reflection coefficients into            */
/*   predictor coefficients.                         */
/*****************************************************/
void StepUp(int fOrder,double* kArray,double* aArray) {

	double aTemp[TNS_MAX_ORDER+2];
	int i,order;

	aArray[0]=1.0;
	aTemp[0]=1.0;
	for (order=1;order<=fOrder;order++) {
		aArray[order]=0.0;
		for (i=1;i<=order;i++) {
			aTemp[i] = aArray[i] + kArray[order]*aArray[order-i];
		}
		for (i=1;i<=order;i++) {
			aArray[i]=aTemp[i];
		}
	}
}


/* Former test code....
main() {
	int i;
	double k[TNS_MAX_ORDER+1];
	double a[TNS_MAX_ORDER+1];
	int ind[TNS_MAX_ORDER+1];
	double e;
	#define DATA_SIZE 1024
	double testData[DATA_SIZE];

	for (i=0;i<DATA_SIZE;i++) {
		testData[i]=sin((double)i/10*2*PI);
		testData[i]+=sin((double)i/5*2*PI);
	}

	e=LevinsonDurbin(TNS_MAX_ORDER,DATA_SIZE,testData,k);
	for (i=0;i<=TNS_MAX_ORDER;i++) {
		fprintf(stdout,"k[%d]=%f\n",i,k[i]);
	}
	fprintf(stdout,"Filter order is:%d\n",TruncateCoeffs(TNS_MAX_ORDER,DEF_TNS_COEFF_THRESH,k));
	for (i=0;i<=TNS_MAX_ORDER;i++) {
		fprintf(stdout,"k[%d]=%f\n",i,k[i]);
	}
	QuantizeReflectionCoeffs(TNS_MAX_ORDER,16,k,ind);
    StepUp(TNS_MAX_ORDER,k,a);

	return 0;
}

*/