filterkit.c

音频信号的重采样程序,如44.1K的WAV转换成采样频率为48K的WAV.

    
    fscanf(fp, "ScaleFactor ");
    if (1 != fscanf(fp,"%d",&temp))
      return(2);
    *LpScl = temp;
    
    fscanf(fp, "\nLength ");
    if (1 != fscanf(fp,"%d",&temp))
      return(3);
    *Nwing = temp;
    
    fscanf(fp, "\nNmult ");
    if (1 != fscanf(fp,"%d",&temp))
      return(4);
    *Nmult = temp;
    
    Imp = (HWORD *) malloc(*Nwing * sizeof(HWORD));
    
    fscanf(fp, "\nCoeffs:\n");
    for (i=0; i<*Nwing; i++)  { /* Get array of 16-bit filter coefficients */
	fscanf(fp, "%d\n", &temp);
	Imp[i] = temp;
    }    

    ImpD = (HWORD *) malloc(*Nwing * sizeof(HWORD));
    
    fscanf(fp, "\nDifferences:\n");
    for (i=0; i<*Nwing; i++)  { /* Get array of 16bit filter coeff differences */
	fscanf(fp, "%d\n", &temp);
	ImpD[i] = temp;
    }    
    
    fclose(fp);
    if (!filterFile || !(*filterFile))
      free(fname);
    *ImpP = Imp;
    *ImpDP = ImpD;
    return(0);
}


WORD FilterUp(HWORD Imp[], HWORD ImpD[], 
		     UHWORD Nwing, BOOL Interp,
		     HWORD *Xp, HWORD Ph, HWORD Inc)
{
    HWORD *Hp, *Hdp = NULL, *End;
    HWORD a = 0;
    WORD v, t;
    
    v=0;
    Hp = &Imp[Ph>>Na];
    End = &Imp[Nwing];
    if (Interp) {
	Hdp = &ImpD[Ph>>Na];
	a = Ph & Amask;
    }
    if (Inc == 1)		/* If doing right wing...              */
    {				/* ...drop extra coeff, so when Ph is  */
	End--;			/*    0.5, we don't do too many mult's */
	if (Ph == 0)		/* If the phase is zero...           */
	{			/* ...then we've already skipped the */
	    Hp += Npc;		/*    first sample, so we must also  */
	    Hdp += Npc;		/*    skip ahead in Imp[] and ImpD[] */
	}
    }
    if (Interp)
      while (Hp < End) {
	  t = *Hp;		/* Get filter coeff */
	  t += (((WORD)*Hdp)*a)>>Na; /* t is now interp'd filter coeff */
	  Hdp += Npc;		/* Filter coeff differences step */
	  t *= *Xp;		/* Mult coeff by input sample */
	  if (t & (1<<(Nhxn-1)))  /* Round, if needed */
	    t += (1<<(Nhxn-1));
	  t >>= Nhxn;		/* Leave some guard bits, but come back some */
	  v += t;			/* The filter output */
	  Hp += Npc;		/* Filter coeff step */
	  Xp += Inc;		/* Input signal step. NO CHECK ON BOUNDS */
      } 
    else 
      while (Hp < End) {
	  t = *Hp;		/* Get filter coeff */
	  t *= *Xp;		/* Mult coeff by input sample */
	  if (t & (1<<(Nhxn-1)))  /* Round, if needed */
	    t += (1<<(Nhxn-1));
	  t >>= Nhxn;		/* Leave some guard bits, but come back some */
	  v += t;			/* The filter output */
	  Hp += Npc;		/* Filter coeff step */
	  Xp += Inc;		/* Input signal step. NO CHECK ON BOUNDS */
      }
    return(v);
}

WORD FilterUD( HWORD Imp[], HWORD ImpD[],
		     UHWORD Nwing, BOOL Interp,
		     HWORD *Xp, HWORD Ph, HWORD Inc, UHWORD dhb)
{
    HWORD a;
    HWORD *Hp, *Hdp, *End;
    WORD v, t;
    UWORD Ho;
    
    v=0;
    Ho = (Ph*(UWORD)dhb)>>Np;
    End = &Imp[Nwing];
    if (Inc == 1)		/* If doing right wing...              */
    {				/* ...drop extra coeff, so when Ph is  */
	End--;			/*    0.5, we don't do too many mult's */
	if (Ph == 0)		/* If the phase is zero...           */
	  Ho += dhb;		/* ...then we've already skipped the */
    }				/*    first sample, so we must also  */
				/*    skip ahead in Imp[] and ImpD[] */
    if (Interp)
      while ((Hp = &Imp[Ho>>Na]) < End) {
	  t = *Hp;		/* Get IR sample */
	  Hdp = &ImpD[Ho>>Na];  /* get interp (lower Na) bits from diff table*/
	  a = Ho & Amask;	/* a is logically between 0 and 1 */
	  t += (((WORD)*Hdp)*a)>>Na; /* t is now interp'd filter coeff */
	  t *= *Xp;		/* Mult coeff by input sample */
	  if (t & 1<<(Nhxn-1))	/* Round, if needed */
	    t += 1<<(Nhxn-1);
	  t >>= Nhxn;		/* Leave some guard bits, but come back some */
	  v += t;			/* The filter output */
	  Ho += dhb;		/* IR step */
	  Xp += Inc;		/* Input signal step. NO CHECK ON BOUNDS */
      }
    else 
      while ((Hp = &Imp[Ho>>Na]) < End) {
	  t = *Hp;		/* Get IR sample */
	  t *= *Xp;		/* Mult coeff by input sample */
	  if (t & 1<<(Nhxn-1))	/* Round, if needed */
	    t += 1<<(Nhxn-1);
	  t >>= Nhxn;		/* Leave some guard bits, but come back some */
	  v += t;			/* The filter output */
	  Ho += dhb;		/* IR step */
	  Xp += Inc;		/* Input signal step. NO CHECK ON BOUNDS */
      }
    return(v);
}

/*
 * double zerox(Data, Factor)
 * HWORD *Data;
 * double Factor;
 *    Given a pointer into a sound sample, this function uses a low-pass
 * filter to estimate the x coordinate of the zero-crossing which must ocurr
 * between Data[0] and Data[1].  This value is returned as the value of the
 * function.  A return value of -100 indicates there was no zero-crossing in
 * the x interval [-1,2].  Factor is the resampling factor: Rate(out) /
 * Rate(in).  Nmult (which determines which filter is used) is passed the
 * zerox's initialization routine: initZerox(Nmult)
 *                                 UHWORD Nmult;
 */

static UHWORD LpScl, Nmult, Nwing;
static HWORD *Imp;
static HWORD *ImpD;

/* ERROR return values:
 *   0 - no error
 *   1 - Nmult is even (should be odd)
 *   2 - filter file not found
 *   3 - invalid ScaleFactor in input file
 *   4 - invalid Length in file
 *   5 - invalid Nmult in file
 */
int initZerox(UHWORD tempNmult)
{
   int err;

   /* Check for illegal input values */
   if (!(tempNmult % 2))
      return(1);
   err = readFilter(NULL, (HWORD **)&Imp, (HWORD **)&ImpD, 
			&LpScl, &tempNmult, &Nwing);
   if (err)
      return(1+err);

   Nmult = tempNmult;
   return(0);
}

#define MAXITER 64
#define ZeroxEPSILON (1E-4)
#define ZeroxMAXERROR (5.0)

double zerox(HWORD *Data, double Factor)
{
   double x, out;
   double lo, hi;
   double dh;
   UWORD dhb;
   WORD v;
   int i;

   if (!Data[0])
      return (0.0);
   if (!Data[1])
      return (1.0);

   if (Data[0] < Data[1])
      {
      lo = -1.0;
      hi =  2.0;
      }
   else
      {
      lo =  2.0;
      hi = -1.0;
      }
   dh = (Factor<1) ? (Factor*Npc) : (Npc);
   dhb = dh * (1<<Na) + 0.5;

   for (i=0; i<MAXITER; i++)
      {
      x = (hi+lo)/2.0;
      v  = FilterUD(Imp,ImpD,Nwing,TRUE,Data,  (HWORD)(x*Pmask),    -1,dhb);
      v += FilterUD(Imp,ImpD,Nwing,TRUE,Data+1,(HWORD)((1-x)*Pmask), 1,dhb);
      v >>= Nhg;
      v *= LpScl;
      out = (double)v / (double)(1<<NLpScl);
      if (out < 0.0)
         lo = x;
      else
         hi = x;
      if (ABS(out) <= ZeroxEPSILON)
         return(x);
      }
   printf("|ZeroX Error| x:%g, \t Data[x]:%d, \t Data[x+1]:%d\n",
      x, *Data, *(Data+1));
   printf("|\tABS(out):%g \t EPSILON:%g\n", ABS(out),ZeroxEPSILON);
   if (ABS(out) <= ZeroxMAXERROR)
      return(x);
   return(-100.0);
}


BOOL Query(char *prompt, BOOL deflt, char *help)
{
   char s[80];

   while (TRUE)
      {
      sprintf(s,"\n%s%s", prompt, (*help) ? " (Type ? for help)" : "");
      getstr(s,(deflt)?"yes":"no",s);
      if (*s=='?' && *help)
         printf(help);
      if (*s=='Y' || *s=='y')
         return(TRUE);
      if (*s=='N' || *s=='n')
         return(FALSE);
      }
}


char *GetString(char *prompt, char *deflt, char *help)
{
   static char s[200];

   while (TRUE)
      {
      sprintf(s,"\n%s%s",prompt, (*help) ? " (Type ? for Help)" : "");
      getstr(s,deflt,s);
      if (*s=='?' && *help)
         printf(help);
      else
         return(s);
      }
}


double GetDouble(char *title, double deflt, char *help)
{
   char s[80],sdeflt[80];
   double newval;

   while (TRUE)
      {
      sprintf(s,"\n%s:%s",title, (*help) ? " (Type ? for Help)" : "");
      sprintf(sdeflt,"%g",deflt);
      getstr(s,sdeflt,s);
      if (*s=='?' && *help)
         printf(help);
      else
         {
         if (!sscanf(s,"%lf",&newval))
            return(deflt);
         return(newval);
	 }
      }
}


unsigned short GetUShort(char *title, unsigned short deflt, char *help)
{
   char s[80],sdeflt[80];
   int newval;

   while (TRUE)
      {
      sprintf(s,"\n%s:%s",title, (*help) ? " (Type ? for Help)" : "");
      sprintf(sdeflt,"%d",deflt);
      getstr(s,sdeflt,s);
      if (*s=='?' && *help)
         printf(help);
      else
         {
         if (!sscanf(s,"%d",&newval))
            printf("unchanged (%d)\n",(newval=deflt));
         if (newval < 0)
            printf("Error: value must be >= zero\n");
         else
            return(newval);
	 }
      }
}