filterkit.c

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

/*
 * filterkit.c (library "filterkit.a"):  Kaiser-windowed low-pass filter support.
 */

/* filterkit.c
 *
 * LpFilter() - Calculates the filter coeffs for a Kaiser-windowed low-pass
 *                   filter with a given roll-off frequency.  These coeffs
 *                   are stored into a array of doubles.
 * writeFilter() - Writes a filter to a file.
 * makeFilter() - Calls LpFilter() to create a filter, then scales the double
 *                   coeffs into an array of half words.
 * readFilter() - Reads a filter from a file.
 * FilterUp() - Applies a filter to a given sample when up-converting.
 * FilterUD() - Applies a filter to a given sample when up- or down-
 *                   converting.
 * initZerox() - Initialization routine for the zerox() function.  Must
 *                   be called before zerox() is called.  This routine loads
 *                   the correct filter so zerox() can use it.
 * zerox() - Given a pointer into a sample, finds a zero-crossing on the
 *                   interval [pointer-1:pointer+2] by iteration.
 * Query() - Ask the user for a yes/no question with prompt, default, 
 *                   and optional help.
 * GetUShort() - Ask the user for a unsigned short with prompt, default,
 *                   and optional help.
 * GetDouble() - Ask the user for a double with prompt, default, and
 *                   optional help.
 * GetString() - Ask the user for a string with prompt, default, and
 *                   optional help.
 */

#include "resample.h"
#include "filterkit.h"

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>

/*
 * Getstr() will print the passed "prompt" as a message to the user, and
 * wait for the user to type an input string.  The string is
 * then copied into "answer".  If the user types just a carriage
 * return, then the string "defaultAnswer" is copied into "answer".
 * ???
 * "Answer" and "defaultAnswer" may be the same string, in which case,
 * the defaultAnswer value will be the original contents of "answer".
 * ???
 */
static void getstr(char *prompt, char *defaultAnswer, char *answer)
{
  char *p,s[200];

  printf("%s (defaultAnswer = %s):",prompt,defaultAnswer);
  strcpy(s, prompt);
  if (!(*s))
    strcpy(s, "input:");
  p = s;
  while (*p && *p != '(')
    p++;
  p--;
  while (*p == ' ')
    p--;
  *(++p) = '\0';

  /* gets(answer); */
  fgets(answer,sizeof(answer),stdin);
  answer[strlen(answer)-1] = '\0';

  if (answer[0] == '\0') {
    strcpy(answer, defaultAnswer);
    printf("\t%s = %s\n",s,answer);
  } else {
    printf("\t%s set to %s\n",s,answer);
  }
}


/* LpFilter()
 *
 * reference: "Digital Filters, 2nd edition"
 *            R.W. Hamming, pp. 178-179
 *
 * Izero() computes the 0th order modified bessel function of the first kind.
 *    (Needed to compute Kaiser window).
 *
 * LpFilter() computes the coeffs of a Kaiser-windowed low pass filter with
 *    the following characteristics:
 *
 *       c[]  = array in which to store computed coeffs
 *       frq  = roll-off frequency of filter
 *       N    = Half the window length in number of coeffs
 *       Beta = parameter of Kaiser window
 *       Num  = number of coeffs before 1/frq
 *
 * Beta trades the rejection of the lowpass filter against the transition
 *    width from passband to stopband.  Larger Beta means a slower
 *    transition and greater stopband rejection.  See Rabiner and Gold
 *    (Theory and Application of DSP) under Kaiser windows for more about
 *    Beta.  The following table from Rabiner and Gold gives some feel
 *    for the effect of Beta:
 *
 * All ripples in dB, width of transition band = D*N where N = window length
 *
 *               BETA    D       PB RIP   SB RIP
 *               2.120   1.50  +-0.27      -30
 *               3.384   2.23    0.0864    -40
 *               4.538   2.93    0.0274    -50
 *               5.658   3.62    0.00868   -60
 *               6.764   4.32    0.00275   -70
 *               7.865   5.0     0.000868  -80
 *               8.960   5.7     0.000275  -90
 *               10.056  6.4     0.000087  -100
 */


#define IzeroEPSILON 1E-21               /* Max error acceptable in Izero */

static double Izero(double x)
{
   double sum, u, halfx, temp;
   int n;

   sum = u = n = 1;
   halfx = x/2.0;
   do {
      temp = halfx/(double)n;
      n += 1;
      temp *= temp;
      u *= temp;
      sum += u;
      } while (u >= IzeroEPSILON*sum);
   return(sum);
}


void LpFilter(double c[], int N, double frq, double Beta, int Num)
{
   double IBeta, temp, inm1;
   int i;

   /* Calculate ideal lowpass filter impulse response coefficients: */
   c[0] = 2.0*frq;
   for (i=1; i<N; i++) {
       temp = PI*(double)i/(double)Num;
       c[i] = sin(2.0*temp*frq)/temp; /* Analog sinc function, cutoff = frq */
   }

   /* 
    * Calculate and Apply Kaiser window to ideal lowpass filter.
    * Note: last window value is IBeta which is NOT zero.
    * You're supposed to really truncate the window here, not ramp
    * it to zero. This helps reduce the first sidelobe. 
    */
   IBeta = 1.0/Izero(Beta);
   inm1 = 1.0/((double)(N-1));
   for (i=1; i<N; i++) {
       temp = (double)i * inm1;
       c[i] *= Izero(Beta*sqrt(1.0-temp*temp)) * IBeta;
   }
}


/* Write a filter to a file
 *    Filter file format:
 *       file name: "F" Nmult "T" Nhc ".filter"
 *       1st line:  the string "ScaleFactor" followed by its value.
 *       2nd line:  the string "Length" followed by Nwing's value.
 *       3rd line:  the string "Nmult" followed by Nmult's value.
 *       4th line:  the string "Coeffs:" on a separate line.
 *       following lines:  Nwing number of 16-bit impulse response values
 *          in the right wing of the impulse response (the Imp[] array).
 *         (Nwing is equal to Npc*(Nmult+1)/2+1, where Npc is defined in the
 *         file "resample.h".)  Each coefficient is on a separate line.
 *       next line:  the string "Differences:" on a separate line.
 *       following lines:  Nwing number of 16-bit impulse-response
 *          successive differences:  ImpD[i] = Imp[i+1] - Imp[i].
 * ERROR codes:
 *   0 - no error
 *   1 - could not open file
 */

int writeFilter(HWORD Imp[], HWORD ImpD[], UHWORD LpScl, UHWORD Nmult, 
		UHWORD Nwing)
{
   char fname[30];
   FILE *fp;
   int i;

   sprintf(fname, "F%dT%d.filter", Nmult, Nhc);
   fp = fopen(fname, "w");
   if (!fp)
      return(1);
   fprintf(fp, "ScaleFactor %d\n", LpScl);
   fprintf(fp, "Length %d\n", Nwing);
   fprintf(fp, "Nmult %d\n", Nmult);
   fprintf(fp, "Coeffs:\n");
   for (i=0; i<Nwing; i++)   /* Put array of 16-bit filter coefficients */
      fprintf(fp, "%d\n", Imp[i]);
   fprintf(fp, "Differences:\n");
   for (i=0; i<Nwing; i++)   /* Put array of 16-bit filter coeff differences */
      fprintf(fp, "%d\n", ImpD[i]);
   fclose(fp);
   printf("Wrote filter file '%s' in current directory.\n",fname);
   return(0);
}


/* ERROR return codes:
 *    0 - no error
 *    1 - Nwing too large (Nwing is > MAXNWING)
 *    2 - Froll is not in interval [0:1)
 *    3 - Beta is < 1.0
 *    4 - LpScl will not fit in 16-bits
 *
 * Made following global to avoid stack problems in Sun3 compilation: */

#define MAXNWING   8192
static double ImpR[MAXNWING];

int makeFilter(HWORD Imp[], HWORD ImpD[], UHWORD *LpScl, UHWORD Nwing,
	       double Froll, double Beta)
{
   double DCgain, Scl, Maxh;
   HWORD Dh;
   int i, temp;

   if (Nwing > MAXNWING)                      /* Check for valid parameters */
      return(1);
   if ((Froll<=0) || (Froll>1))
      return(2);
   if (Beta < 1)
      return(3);

   /* 
    * Design Kaiser-windowed sinc-function low-pass filter 
    */
   LpFilter(ImpR, (int)Nwing, 0.5*Froll, Beta, Npc); 

   /* Compute the DC gain of the lowpass filter, and its maximum coefficient
    * magnitude. Scale the coefficients so that the maximum coeffiecient just
    * fits in Nh-bit fixed-point, and compute LpScl as the NLpScl-bit (signed)
    * scale factor which when multiplied by the output of the lowpass filter
    * gives unity gain. */
   DCgain = 0;
   Dh = Npc;                       /* Filter sampling period for factors>=1 */
   for (i=Dh; i<Nwing; i+=Dh)
      DCgain += ImpR[i];
   DCgain = 2*DCgain + ImpR[0];    /* DC gain of real coefficients */

   for (Maxh=i=0; i<Nwing; i++)
      Maxh = MAX(Maxh, fabs(ImpR[i]));

   Scl = ((1<<(Nh-1))-1)/Maxh;     /* Map largest coeff to 16-bit maximum */
   temp = fabs((1<<(NLpScl+Nh))/(DCgain*Scl));
   if (temp >= 1<<16)
      return(4);                   /* Filter scale factor overflows UHWORD */
   *LpScl = temp;

   /* Scale filter coefficients for Nh bits and convert to integer */
   if (ImpR[0] < 0)                /* Need pos 1st value for LpScl storage */
      Scl = -Scl;
   for (i=0; i<Nwing; i++)         /* Scale them */
      ImpR[i] *= Scl;
   for (i=0; i<Nwing; i++)         /* Round them */
      Imp[i] = ImpR[i] + 0.5;

   /* ImpD makes linear interpolation of the filter coefficients faster */
   for (i=0; i<Nwing-1; i++)
      ImpD[i] = Imp[i+1] - Imp[i];
   ImpD[Nwing-1] = - Imp[Nwing-1];      /* Last coeff. not interpolated */

   return(0);
}


/* Read-in a filter
 *    Filter file format:
 *       file name: "F" Nmult "T" Nhc ".filter"
 *       1st line:  the string "ScaleFactor" followed by its value.
 *       2nd line:  the string "Length" followed by Nwing's value.
 *       3rd line:  the string "Coeffs:" on separate line.
 *       Nwing number of 16-bit impulse response values in the right
 *          wing of the impulse response.  (Length=Npc*(Nmult+1)/2+1,
 *          where originally Npc=2^9, and Nmult=13.)   Each on separate line.
 *       The string "Differences:" on separate line.
 *       Nwing number of 16-bit impulse-response successive differences:
 *          ImpDiff[i] = Imp[i+1] - Imp[i].
 *
 * ERROR return codes:
 *    0 - no error
 *    1 - file not found
 *    2 - invalid ScaleFactor in file
 *    3 - invalid Length in file
 *    4 - invalid Nmult in file
 */
int readFilter(char *filterFile, HWORD **ImpP, HWORD **ImpDP, UHWORD *LpScl, 
	       UHWORD *Nmult, UHWORD *Nwing)
{
    char *fname;
    FILE *fp;
    int i, temp;
    HWORD *Imp,*ImpD;
    
    if (!filterFile || !(*filterFile)) {
	fname = (char *) malloc(32);
	if ((*Nmult)>0 && ((*Nmult)&1))
	  sprintf(fname, "F%dT%d.filter", *Nmult, Nhc);
	else
	  sprintf(fname, "F65dT%d.filter", Nhc);
    } else
      fname = filterFile;
    
    fp = fopen(fname, "r");
    if (fp == NULL)
      return(1);