#!/bin/sh
# This is part 02 of resample
# ============= filterkit.c ==============
if test -f 'filterkit.c' -a X"$1" != X"-c"; then
	echo 'x - skipping filterkit.c (File already exists)'
else
echo 'x - extracting filterkit.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'filterkit.c' &&
X
X/*
X * FILE: filterkit.c (source for library "filterkit.a")
X *   BY: Christopher Lee Fraley (cf0v@andrew.cmu.edu)
X * DESC: Makes a Kaiser-windowed low-pass filter.
X * DATE: 7-JUN-88
X */
X
Xchar filterkitVERSION[] = "2.0  (7-JUN-88  3:30pm)";
X
X/* filterkit.c
X *
X *      The function makeFilter(), FilterUp(), and FilterUD(), were originally
X * written by Julius O. Smith in SAIL, and were translated into C.
X *
X * LpFilter() - Calculates the filter coeffs for a Kaiser-windowed low-pass
X *                   filter with a given roll-off frequency.  These coeffs
X *                   are stored into a array of doubles.
X * writeFilter() - Writes a filter to a file.
X * makeFilter() - Calls LpFilter() to create a filter, then scales the double
X *                   coeffs into an array of half words.
X * readFilter() - Reads a filter from a file.
X * FilterUp() - Applies a filter to a given sample when up-converting.
X * FilterUD() - Applies a filter to a given sample when up- or down-
X *                   converting.
X * initZerox() - Initialization routine for the zerox() function.  Must
X *                   be called before zerox() is called.  This routine loads
X *                   the correct filter so zerox() can use it.
X * zerox() - Given a pointer into a sample, finds a zero-crossing on the
X *                   interval [pointer-1:pointer+2] by iteration.
X * Query() - Ask the user for a yes/no question with prompt, default, 
X *                   and optional help.
X * GetUShort() - Ask the user for a unsigned short with prompt, default,
X *                   and optional help.
X * GetDouble() - Ask the user for a double with prompt, default, and
X *                   optional help.
X * GetString() - Ask the user for a string with prompt, default, and
X *                   optional help.
X */
X
X
X#include <stdio.h>
X#include <math.h>
X#include <string.h>
X#include "stdefs.h"
X#include "resample.h"
X#include "filterkit.h"
X#include "protos.h"
X
X
X
X/* LpFilter()
X *
X * reference: "Digital Filters, 2nd edition"
X *            R.W. Hamming, pp. 178-179
X *
X * Izero() computes the 0th order modified bessel function of the first kind.
X *    (Needed to compute Kaiser window).
X *
X * LpFilter() computes the coeffs of a Kaiser-windowed low pass filter with
X *    the following characteristics:
X *
X *       c[]  = array in which to store computed coeffs
X *       frq  = roll-off frequency of filter
X *       N    = Half the window length in number of coeffs
X *       Beta = parameter of Kaiser window
X *       Num  = number of coeffs before 1/frq
X *
X * Beta trades the rejection of the lowpass filter against the transition
X *    width from passband to stopband.  Larger Beta means a slower
X *    transition and greater stopband rejection.  See Rabiner and Gold
X *    (Theory and Application of DSP) under Kaiser windows for more about
X *    Beta.  The following table from Rabiner and Gold gives some feel
X *    for the effect of Beta:
X *
X * All ripples in dB, width of transition band = D*N where N = window length
X *
X *               BETA    D       PB RIP   SB RIP
X *               2.120   1.50  +-0.27      -30
X *               3.384   2.23    0.0864    -40
X *               4.538   2.93    0.0274    -50
X *               5.658   3.62    0.00868   -60
X *               6.764   4.32    0.00275   -70
X *               7.865   5.0     0.000868  -80
X *               8.960   5.7     0.000275  -90
X *               10.056  6.4     0.000087  -100
X */
X
X
X
X#define IzeroEPSILON 1E-21               /* Max error acceptable in Izero */
X
Xdouble Izero(x)
Xdouble x;
X{
X   double sum, u, halfx, temp;
X   int n;
X
X   sum = u = n = 1;
X   halfx = x/2.0;
X   do {
X      temp = halfx/(double)n;
X      n += 1;
X      temp *= temp;
X      u *= temp;
X      sum += u;
X      } while (u >= IzeroEPSILON*sum);
X   return(sum);
X}
X
X
Xvoid LpFilter(c,N,frq,Beta,Num)
Xdouble c[], frq, Beta;
Xint N, Num;
X{
X   double IBeta, temp;
X   int i;
X
X   /* Calculate filter coeffs: */
X   c[0] = 2.0*frq;
X   for (i=1; i<N; i++)
X      {
X      temp = PI*(double)i/(double)Num;
X      c[i] = sin(2.0*temp*frq)/temp;
X      }
X
X   /* Calculate and Apply Kaiser window to filter coeffs: */
X   IBeta = 1.0/Izero(Beta);
X   for (i=1; i<N; i++)
X      {
X      temp = (double)i / ((double)N * (double)1.0);
X      c[i] *= Izero(Beta*sqrt(1.0-temp*temp)) * IBeta;
X      }
X}
X
X
X
X/* Write a filter to a file
X *    Filter file format:
X *       file name: "F" Nmult "T" Nhc ".filter"
X *       1st line:  the string "ScaleFactor" followed by its value.
X *       2nd line:  the string "Length" followed by Nwing's value.
X *       3rd line:  the string "Coeffs:" on a separate line.
X *       following lines:  Nwing number of 16-bit impulse response values
X *          in the right wing of the impulse response (the Imp[] array).
X *         (Nwing is equal to Npc*(Nmult+1)/2+1, where Npc is defined in the
X *         file "resample.h".)  Each coefficient is on a separate line.
X *       next line:  the string "Differences:" on a separate line.
X *       following lines:  Nwing number of 16-bit impulse-response
X *          successive differences:  ImpD[i] = Imp[i+1] - Imp[i].
X * ERROR codes:
X *   0 - no error
X *   1 - could not open file
X */
X
Xint writeFilter(Imp, ImpD, LpScl, Nmult, Nwing)
XHWORD Imp[], ImpD[];
XUHWORD LpScl, Nmult, Nwing;
X{
X   char fname[30];
X   FILE *fp;
X   int i;
X
X   sprintf(fname, "F%dT%d.filter", Nmult, Nhc);
X   fp = fopen(fname, "w");
X   if (!fp)
X      return(1);
X   fprintf(fp, "ScaleFactor %d\n", LpScl);
X   fprintf(fp, "Length %d\n", Nwing);
X   fprintf(fp, "Coeffs:\n");
X   for (i=0; i<Nwing; i++)   /* Put array of 16-bit filter coefficients */
X      fprintf(fp, "%d\n", Imp[i]);
X   fprintf(fp, "Differences:\n");
X   for (i=0; i<Nwing; i++)   /* Put array of 16-bit filter coeff differences */
X      fprintf(fp, "%d\n", ImpD[i]);
X   fclose(fp);
X   return(0);
X}
X
X
X
X/* ERROR return codes:
X *    0 - no error
X *    1 - Nwing too large (Nwing is > MAXNWING)
X *    2 - Froll is not in interval [0:1)
X *    3 - Beta is < 1.0
X *    4 - LpScl will not fit in 16-bits
X *
X * Made following global to avoid stack problems in Sun3 compilation: */
Xstatic double ImpR[MAXNWING];
X
Xint makeFilter(Imp, ImpD, LpScl, Nwing, Froll, Beta)
XHWORD Imp[], ImpD[];
XUHWORD *LpScl, Nwing;
Xdouble Froll, Beta;
X{
X   double DCgain, Scl, Maxh;
X   HWORD Dh;
X   int i, temp;
X
X   if (Nwing > MAXNWING)                      /* Check for valid parameters */
X      return(1);
X   if ((Froll<=0) || (Froll>1))
X      return(2);
X   if (Beta < 1)
X      return(3);
X
X   LpFilter(ImpR, (int)Nwing, Froll, Beta, Npc); /* Design a Kaiser-window */
X                                                 /* Sinc low-pass filter */
X
X   /* Compute the DC gain of the lowpass filter, and its maximum coefficient
X    * magnitude. Scale the coefficients so that the maximum coeffiecient just
X    * fits in Nh-bit fixed-point, and compute LpScl as the NLpScl-bit (signed)
X    * scale factor which when multiplied by the output of the lowpass filter
X    * gives unity gain. */
X   DCgain = 0;
X   Dh = Npc;                       /* Filter sampling period for factors>=1 */
X   for (i=Dh; i<Nwing; i+=Dh)
X      DCgain += ImpR[i];
X   DCgain = 2*DCgain + ImpR[0];    /* DC gain of real coefficients */
X
X   for (Maxh=i=0; i<Nwing; i++)
X      Maxh = MAX(Maxh, fabs(ImpR[i]));
X
X   Scl = ((1<<(Nh-1))-1)/Maxh;     /* Map largest coeff to 16-bit maximum */
X   temp = fabs((1<<(NLpScl+Nh))/(DCgain*Scl));
X   if (temp >= 1<<16)
X      return(4);                   /* Filter scale factor overflows UHWORD */
X   *LpScl = temp;
X
X   /* Scale filter coefficients for Nh bits and convert to integer */
X   if (ImpR[0] < 0)                /* Need pos 1st value for LpScl storage */
X      Scl = -Scl;
X   for (i=0; i<Nwing; i++)         /* Scale them */
X      ImpR[i] *= Scl;
X   for (i=0; i<Nwing; i++)         /* Round them */
X      Imp[i] = ImpR[i] + 0.5;
X
X   /* ImpD makes linear interpolation of the filter coefficients faster */
X   for (i=0; i<Nwing-1; i++)
X      ImpD[i] = Imp[i+1] - Imp[i];
X   ImpD[Nwing-1] = - Imp[Nwing-1];      /* Last coeff. not interpolated */
X
X   return(0);
X}
X
X
X
X/* Read-in a filter
X *    Filter file format:
X *       file name: "F" Nmult "T" Nhc ".filter"
X *       1st line:  the string "ScaleFactor" followed by its value.
X *       2nd line:  the string "Length" followed by Nwing's value.
X *       3rd line:  the string "Coeffs:" on separate line.
X *       Nwing number of 16-bit impulse response values in the right
X *          wing of the impulse response.  (Length=Npc*(Nmult+1)/2+1,
X *          where originally Npc=2^9, and Nmult=13.)   Each on separate line.
X *       The string "Differences:" on separate line.
X *       Nwing number of 16-bit impulse-response successive differences:
X *          ImpDiff[i] = Imp[i+1] - Imp[i].
X *
X * ERROR return codes:
X *    0 - no error
X *    1 - file not found
X *    2 - invalid ScaleFactor in file
X *    3 - invalid Length in file
X */
Xint readFilter(Imp, ImpD, LpScl, Nmult, Nwing)
XHWORD Imp[], ImpD[];
XUHWORD *LpScl, Nmult, *Nwing;
X{
X   char fname[30];
X   FILE *fp;
X   int i, temp;
X
X   sprintf(fname, "F%dT%d.filter", Nmult, Nhc);
X   fp = fopen(fname, "r");
X   if (fp == NULL)
X      return(1);
X
X   fscanf(fp, "ScaleFactor ");
X   if (1 != fscanf(fp,"%d",&temp))
X      return(2);
X   *LpScl = temp;
X
X   fscanf(fp, "\nLength ");
X   if (1 != fscanf(fp,"%d",&temp))
X      return(3);
X   *Nwing = temp;
X
X   fscanf(fp, "\nCoeffs:\n");
X   for (i=0; i<*Nwing; i++)  /* Get array of 16-bit filter coefficients */
X      fscanf(fp, "%d\n", &Imp[i]);
X
X   fscanf(fp, "\nDifferences:\n");
X   for (i=0; i<*Nwing; i++)  /* Get array of 16-bit filter coeff differences */
X      fscanf(fp, "%d\n", &ImpD[i]);
X
X   fclose(fp);
X   return(0);
X}
X
X
X
X
XWORD FilterUp(Imp, ImpD, Nwing, Interp, Xp, Ph, Inc)
XHWORD Imp[], ImpD[], *Xp, Inc, Ph;
XUHWORD Nwing;
XBOOL Interp;
X{
X   HWORD a, *Hp, *Hdp, *End;
X   WORD v, t;
X
X   v=0;
X   Hp = &Imp[Ph>>Na];
X   End = &Imp[Nwing];
X   if (Interp)
X      {
X      Hdp = &ImpD[Ph>>Na];
X      a = Ph & Amask;
X      }
X   if (Inc == 1)                     /* If doing right wing...              */
X      {                              /* ...drop extra coeff, so when Ph is  */
X      End--;                         /*    0.5, we don't do too many mult's */
X      if (Ph == 0)                   /* If the phase is zero...           */
X         {                           /* ...then we've already skipped the */
X         Hp += Npc;                  /*    first sample, so we must also  */
X         Hdp += Npc;                 /*    skip ahead in Imp[] and ImpD[] */
X         }
X      }
X   while (Hp < End)
X      {
X      t = *Hp;                       /* Get filter coeff */
X      if (Interp)
X         {
X         t += (((WORD)*Hdp)*a)>>Na;  /* t is now interp'd filter coeff */
X         Hdp += Npc;                 /* Filter coeff differences step */
X	 }
X      t *= *Xp;      /* Mult coeff by input sample */
X      if (t & 1<<(Nhxn-1))  /* Round, if needed */
X         t += 1<<(Nhxn-1);
X      t >>= Nhxn;    /* Leave some guard bits, but come back some */
X      v += t;        /* The filter output */
X      Hp += Npc;     /* Filter coeff step */
X      Xp += Inc;     /* Input signal step. NO CHECK ON ARRAY BOUNDS */
X      }
X   return(v);
X}
X
X
X
X
XWORD FilterUD(Imp, ImpD, Nwing, Interp, Xp, Ph, Inc, dhb)
XHWORD Imp[], ImpD[], *Xp, Ph, Inc;
XUHWORD Nwing, dhb;
XBOOL Interp;
X{
X   HWORD a, *Hp, *Hdp, *End;
X   WORD v, t;
X   UWORD Ho;
X
X   v=0;
X   Ho = (Ph*(UWORD)dhb)>>Np;
X   End = &Imp[Nwing];
X   if (Inc == 1)                     /* If doing right wing...              */
X      {                              /* ...drop extra coeff, so when Ph is  */
X      End--;                         /*    0.5, we don't do too many mult's */
X      if (Ph == 0)                   /* If the phase is zero...           */
X         Ho += dhb;                  /* ...then we've already skipped the */
X      }                              /*    first sample, so we must also  */
X                                     /*    skip ahead in Imp[] and ImpD[] */
X   while ((Hp = &Imp[Ho>>Na]) < End)
X      {
X      t = *Hp;       /* Get IR sample */
X      if (Interp)
X         {
X         Hdp = &ImpD[Ho>>Na]; /* get interp (lower Na) bits from diff table */
X         a = Ho & Amask;                  /* a is logically between 0 and 1 */
X         t += (((WORD)*Hdp)*a)>>Na;       /* t is now interp'd filter coeff */
X	 }
X      t *= *Xp;      /* Mult coeff by input sample */
X      if (t & 1<<(Nhxn-1))  /* Round, if needed */
X         t += 1<<(Nhxn-1);
X      t >>= Nhxn;    /* Leave some guard bits, but come back some */
X      v += t;        /* The filter output */
X      Ho += dhb;     /* IR step */
X      Xp += Inc;     /* Input signal step. NO CHECK ON ARRAY BOUNDS */
X      }
X   return(v);
X}
X
X
X
X/*
X * double zerox(Data, Factor)
X * HWORD *Data;
X * double Factor;
X *    Given a pointer into a sound sample, this function uses a low-pass
X * filter to estimate the x coordinate of the zero-crossing which must ocurr
X * between Data[0] and Data[1].  This value is returned as the value of the
X * function.  A return value of -100 indicates there was no zero-crossing in
X * the x interval [-1,2].  Factor is the resampling factor: Rate(out) /
X * Rate(in).  Nmult (which determines which filter is used) is passed the
X * zerox's initialization routine: initZerox(Nmult)
X *                                 UHWORD Nmult;
X */
X
Xstatic UHWORD LpScl, Nmult, Nwing;
Xstatic HWORD Imp[MAXNWING];
Xstatic HWORD ImpD[MAXNWING];
X
X/* ERROR return values:
X *   0 - no error
X *   1 - Nmult is even (should be odd)
X *   2 - filter file not found
X *   3 - invalid ScaleFactor in input file
X *   4 - invalid Length in file
X */
Xint initZerox(tempNmult)
XUHWORD tempNmult;
X{
X   int err;
X
X   /* Check for illegal input values */
X   if (!(tempNmult % 2))
X      return(1);
X   if (err = readFilter(Imp, ImpD, &LpScl, tempNmult, &Nwing))
X      return(1+err);
X
X   Nmult = tempNmult;
X   return(0);
X}
X
X
X
X#define MAXITER 64
X#define ZeroxEPSILON (1E-4)
X#define ZeroxMAXERROR (5.0)
X
Xdouble zerox(Data, Factor)
XHWORD *Data;
Xdouble Factor;
X{
X   double x, out;
X   double lo, hi;
X   double dh;
X   UWORD dhb;
X   WORD v;
X   int i;
X
X   if (!Data[0])
X      return (0.0);
X   if (!Data[1])
X      return (1.0);
X
X   if (Data[0] < Data[1])
X      {
X      lo = -1.0;
X      hi =  2.0;
X      }
X   else
X      {
X      lo =  2.0;
X      hi = -1.0;
X      }
X   dh = (Factor<1) ? (Factor*Npc) : (Npc);
X   dhb = dh * (1<<Na) + 0.5;
X
X   for (i=0; i<MAXITER; i++)
X      {
X      x = (hi+lo)/2.0;
X      v  = FilterUD(Imp,ImpD,Nwing,TRUE,Data,  (HWORD)(x*Pmask),    -1,dhb);
X      v += FilterUD(Imp,ImpD,Nwing,TRUE,Data+1,(HWORD)((1-x)*Pmask), 1,dhb);
X      v >>= Nhg;
X      v *= LpScl;
X      out = (double)v / (double)(1<<NLpScl);
X      if (out < 0.0)
X         lo = x;
X      else
X         hi = x;
X      if (ABS(out) <= ZeroxEPSILON)
X         return(x);
X      }
X   printf("|ZeroX Error| x:%g, \t Data[x]:%d, \t Data[x+1]:%d\n",
X      x, *Data, *(Data+1));
X   printf("|\tABS(out):%g \t EPSILON:%g\n", ABS(out),ZeroxEPSILON);
X   if (ABS(out) <= ZeroxMAXERROR)
X      return(x);
X   return(-100.0);
X}
X
X
X
X
XBOOL Query(prompt, deflt, help)
Xchar *help, *prompt;
XBOOL deflt;
X{
X   char s[80];
X
X   while (TRUE)
X      {
X      sprintf(s,"\n%s%s", prompt, (*help) ? " (Type ? for help)" : "");
X      getstr(s,(deflt)?"yes":"no",s);
X      if (*s=='?' && *help)
X         printf(help);
X      if (*s=='Y' || *s=='y')
X         return(TRUE);
X      if (*s=='N' || *s=='n')
X         return(FALSE);
X      }
X}
X
X
Xchar *GetString(prompt, deflt, help)
Xchar *help, *prompt, *deflt;
X{
X   char s[80];
X
X   while (TRUE)
X      {
X      sprintf(s,"\n%s%s",prompt, (*help) ? " (Type ? for Help)" : "");
X      getstr(s,deflt,s);
X      if (*s=='?' && *help)
X         printf(help);
X      else
X         return(s);
X      }
X}
X
Xgetstr(s1, deflt, s2)
Xchar *s1, *deflt, *s2;
X{
X	write(1, s1, strlen(s1));
X	read(0, s2, 80);
X}
X
X
X
Xdouble GetDouble(title, deflt, help)
Xchar *help, *title;
Xdouble deflt;
X{
X   char s[80],sdeflt[80];
X   double newval;
X
X   while (TRUE)
X      {
X      sprintf(s,"\n%s:%s",title, (*help) ? " (Type ? for Help)" : "");
X      sprintf(sdeflt,"%g",deflt);
X      getstr(s,sdeflt,s);
X      if (*s=='?' && *help)
X         printf(help);
X      else
X         {
X         if (!sscanf(s,"%lf",&newval))
X            return(deflt);
X         return(newval);
X	 }
X      }
X}
X
X
Xunsigned short GetUShort(title, deflt, help)
Xchar *help, *title;
Xunsigned short deflt;
X{
X   char s[80],sdeflt[80];
X   int newval;
X
X   while (TRUE)
X      {
X      sprintf(s,"\n%s:%s",title, (*help) ? " (Type ? for Help)" : "");
X      sprintf(sdeflt,"%d",deflt);
X      getstr(s,sdeflt,s);
X      if (*s=='?' && *help)
X         printf(help);
X      else
X         {
X         if (!sscanf(s,"%d",&newval))
X            printf("unchanged (%d)\n",(newval=deflt));
X         if (newval < 0)
X            printf("Error: value must be >= zero\n");
X         else
X            return(newval);
X	 }
X      }
X}
X
X
SHAR_EOF
chmod 0644 filterkit.c ||
echo 'restore of filterkit.c failed'
Wc_c="`wc -c < 'filterkit.c'`"
test 16885 -eq "$Wc_c" ||
	echo 'filterkit.c: original size 16885, current size' "$Wc_c"
fi
# ============= makefilter.c ==============
if test -f 'makefilter.c' -a X"$1" != X"-c"; then
	echo 'x - skipping makefilter.c (File already exists)'
else
echo 'x - extracting makefilter.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'makefilter.c' &&
X
X/*
X * FILE: makefilter.c
X *   BY: Christopher Lee Fraley (cf0v@spice.cmu.edu or cf0v@andrew.cmu.edu)
X * DESC: Makes a Kaiser-windowed low-pass filter.
X * DATE: 7-JUN-88
X */
X
Xchar makefilterVERSION[] = "2.0  (17-JUN-88  3:00pm)";
X
X#include <stdio.h>
X#include <math.h>
X#include "stdefs.h"
X#include "resample.h"
X#include "filterkit.h"
X
X/* LIBRARIES needed:
X *
X * 1. filterkit
X *       makeFilter()  - designs a Kaiser-windowed low-pass filter
X *       writeFilter() - writes a filter to a standard filter file
X *       GetUShort()   - prompt user for a UHWORD with help
X *       GetDouble()   - prompt user for a double with help
X *
X * 2. math
X */
X
X
X
Xchar NmultHelp[] =
X  "\n   Nmult is the length of the symmetric FIR lowpass filter used\
X   \n   by the sampling rate converter. It must be odd.\
X   \n   This is the number of multiplies per output sample for\
X   \n   up-conversions (Factor>1), and is the number of multiplies\
X   \n   per input sample for down-conversions (Factor<1). Thus if\
X   \n   the rate conversion is Srate2 = Factor*Srate1, then you have\
X   \n   Nmult*Srate1*MAXof(Factor,1) multiplies per second of real time.\
X   \n   Naturally, higher Nmult gives better lowpass-filtering at the\
X   \n   expense of longer compute times. Nmult should be odd because\
X   \n   it is the length of a symmetric FIR filter, and the current\
X   \n   implementation requires a coefficient at the time origin.\n";      
X
Xchar FrollHelp[] =
X  "\n   Froll determines the frequency at which the lowpass filter begins to\
X   \n   roll-off. If Froll=1, then there is no 'guard zone' and the filter\
X   \n   roll-off region will be aliased. If Froll is 0.85, for example, then\
X   \n   the filter begins rolling off at 0.85*Srate/2, so that by Srate/2,\
X   \n   the filter is well down and aliasing is reduced.  Since aliasing\
X   \n   distortion is worse by far than loss of the high-frequency spectral\
X   \n   amplitude, Froll<1 is highly recommended. The default of 0.85\
X   \n   sacrifices the upper 15% of the spectrum as an anti-aliasing guard\
X   \n   zone.\n";
X
Xchar BetaHelp[] =
X  "\n   Beta trades the rejection of the lowpass filter against the\
X   \n   transition width from passband to stopband. Larger Beta means\
X   \n   a slower transition and greater stopband rejection. See Rabiner\
X   \n   and Gold (Th. and App. of DSP) under Kaiser windows for more about\
X   \n   Beta. The following table from Rabiner and Gold gives some feel\
X   \n   for the effect of Beta:\
X   \n\
X   \nAll ripples in dB, width of transition band =D*N, where N= window length\
X   \n\
X   \n               BETA    D       PB RIP   SB RIP\
X   \n               2.120   1.50  +-0.27      -30\
X   \n               3.384   2.23    0.0864    -40\
X   \n               4.538   2.93    0.0274    -50\
X   \n               5.658   3.62    0.00868   -60\
X   \n               6.764   4.32    0.00275   -70\
X   \n               7.865   5.0     0.000868  -80\
X   \n               8.960   5.7     0.000275  -90\
X   \n               10.056  6.4     0.000087  -100\n";
X
X
X
Xmain()
X{
X   HWORD Imp[MAXNWING];               /* Filter coefficients */
X   HWORD ImpD[MAXNWING];              /* ImpD[i] = ImpD[i+1] - ImpD[i] */
X   double Froll, Beta;
X   UHWORD Nmult, Nwing, LpScl;
X   int err;
X
X   printf("\nKaiser-windowed FIR Filter Design\n");
X   printf("Written by:  Chritopher Lee Fraley\n");
X   printf("Version %s\n", makefilterVERSION);
X
X   Nmult = 13;
X   Froll = 0.425;
X   Beta  = 5.7;
X   while (TRUE)
X      {
X      Nmult = GetUHWORD("(Odd) Filter length 'Nmult'", Nmult, NmultHelp);
X      Nwing = Npc*(Nmult+1)/2;     /* # of filter coeffs in right wing */
X      Nwing += Npc/2 + 1;          /* This prevents just missing last coeff */
X                                   /*   for integer conversion factors  */
X      Froll = GetDouble("Normalized Roll-off freq (0<Froll<=1)",
X         Froll, FrollHelp);
X      Beta = GetDouble("Beta", Beta, BetaHelp);
X      printf("\n");
X      if (!(Nmult % 2))
X         printf("Error: Nmult must be odd and greater than zero\n");
X      else if (Nwing > MAXNWING)
X         printf("Error: Nmult too large for current MAXNWING\n");
X      else if ((Froll<=0) || (Froll>1))
X         printf("Error: Roll-off freq must be 0<Froll<=1\n");
X      else if (Beta < 1)
X         printf("Error: Beta must be greater or equal to 1\n");
X      else if (err = makeFilter(Imp, ImpD, &LpScl, Nwing, Froll, Beta))
X         printf("Error: Unable to make filter, err=%d\n", err);
X      else if (err = writeFilter(Imp, ImpD, LpScl, Nmult, Nwing))
X         printf("Error: Unable to write filter, err=%d\n", err);
X      else
X         break;
X      }
X}
X
SHAR_EOF
chmod 0644 makefilter.c ||
echo 'restore of makefilter.c failed'
Wc_c="`wc -c < 'makefilter.c'`"
test 4633 -eq "$Wc_c" ||
	echo 'makefilter.c: original size 4633, current size' "$Wc_c"
fi
# ============= resample.c ==============
if test -f 'resample.c' -a X"$1" != X"-c"; then
	echo 'x - skipping resample.c (File already exists)'
else
echo 'x - extracting resample.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'resample.c' &&
X
X/*
X * FILE: resample.c
X *   BY: Julius Smith (at CCRMA, Stanford U)
X * C BY: translated from SAIL to C by Christopher Lee Fraley
X *          (cf0v@spice.cs.cmu.edu or @andrew.cmu.edu)
X * DATE: 7-JUN-88
X */
X
Xchar resampleVERSION[] = "1.3  (24-JUN-88, 3:00pm)";
X
X/*
X *      The original version of this program in SAIL may be found in:
X *      /../s/usr/mhs/resample  or  /../x/usr/rbd/src/resample
X *
X *      Implement sampling rate conversions by (almost) arbitrary factors.
X *      Based on SRCONV.SAI[SYS,MUS] with new algorithm by JOS&PG.
X *      The program internally uses 16-bit data and 16-bit filter
X *      coefficients.
X *
X *      Reference: "A Flexible Sampling-Rate Conversion Method,"
X *      J. O. Smith and P. Gossett, ICASSP, San Diego, 1984, Pgs 19.4.
X *
X *      * Need overflow detection in Filter() routines.  Also, we want
X *        to saturate instead of wrap-around and report number of clipped
X *        samples.
X */
X
X/* CHANGES from original SAIL program:
X *
X * 1. LpScl is scaled by Factor (when Factor<1) in resample() so this is
X *       done whether the filter was loaded or created.
X * 2. makeFilter() - ImpD[] is created from Imp[] instead of ImpR[], to avoid
X *       problems with round-off errors.
X * 3. makeFilter() - ImpD[Nwing-1] gets NEGATIVE Imp[Nwing-1].
X * 4. SrcU/D() - Switched order of making guard bits (v>>Nhg) and
X *       normalizing.  This was done to prevent overflow.
X */
X
X/* LIBRARIES needed:
X *
X * 1. filterkit
X *       readFilter() - reads standard filter file
X *       FilterUp()   - applies filter to sample when Factor >= 1
X *       FilterUD()   - applies filter to sample for any Factor
X *       Query()      - prompt user for y/n answer with help.
X *       GetDouble()  - prompt user for a double with help.
X *       GetUShort()  - prompt user for a UHWORD with help.
X *
X * 2. math
X */
X
X
X
X
X#include <stdio.h>
X#include <math.h>
X#include <string.h>
X#include "stdefs.h"
X#include "filterkit.h"
X#include "resample.h"
X
X#define IBUFFSIZE 1024                         /* Input buffer size */
X#define OBUFFSIZE (IBUFFSIZE*MAXFACTOR+2)      /* Calc'd out buffer size */
X
X
Xfail(s)
Xchar *s;
X{
X   printf("resample: %s\n\n", s);  /* Display error message  */
X   exit(1);                        /* Exit, indicating error */
X}
X
Xfails(s,s2)
Xchar *s, *s2;
X{
X   printf("resample: ");           /* Display error message  */
X   printf(s,s2);
X   printf("\n\n");
X   exit(1);                        /* Exit, indicating error */
X}
X
X
X
X
X/* Help Declarations */
Xchar FactorHelp[] =
X  "\n   Factor is the amount by which the sampling rate is changed.\
X   \n   If the sampling rate of the input signal is Srate1, then the\
X   \n   sampling rate of the output is Factor*Srate1.\n";
X
Xchar NmultHelp[] =
X  "\n   Nmult is the length of the symmetric FIR lowpass filter used\
X   \n   by the sampling rate converter. It must be odd.\
X   \n   This is the number of multiplies per output sample for\
X   \n   up-conversions (Factor>1), and is the number of multiplies\
X   \n   per input sample for down-conversions (Factor<1). Thus if\
X   \n   the rate conversion is Srate2 = Factor*Srate1, then you have\
X   \n   Nmult*Srate1*MAXof(Factor,1) multiplies per second of real time.\
X   \n   Naturally, higher Nmult gives better lowpass-filtering at the\
X   \n   expense of longer compute times. Nmult should be odd because\
X   \n   it is the length of a symmetric FIR filter, and the current\
X   \n   implementation requires a coefficient at the time origin.\n";      
X
Xchar InterpFiltHelp[] =
X  "\n   By electing to interpolate the filter look-up table,\
X   \n   execution is slowed down but the quality of filtering\
X   \n   is higher. Interpolation results in twice as many\
X   \n   multiply-adds per sample of output.\n";
X
X
X
X/* Global file pointers: */
XFILE *fin, *fout;
X
X
XopenData(argc,argv)
Xint argc;
Xchar *argv[];
X{
X   if (argc != 3)
X      fail("format is 'resample <filein> <fileout>'");
X   if (NULL == (fin = fopen(*++argv,"r")))
X      fails("could not open input file '%s'",*argv);
X   if (NULL == (fout = fopen(*++argv,"w")))
X      fails("could not open output file '%s'",*argv);
X}
X
X
XcloseData()
X{
X   (void) fclose(fin);
X   (void) fclose(fout);
X}
X
X
Xint readData(Data, DataArraySize, Xoff)  /* return: 0 - notDone */
XHWORD Data[];                            /*        >0 - index of last sample */
Xint DataArraySize, Xoff;
X{
X   int Nsamps, Scan;
X   short val;
X   HWORD *DataStart;
X
X   DataStart = Data;
X   Nsamps = DataArraySize - Xoff;   /* Calculate number of samples to get */
X   Data += Xoff;                    /* Start at designated sample number */
X   for (; Nsamps>0; Nsamps--)
X      {
X      Scan = fread(&val, 1, 2, fin);      /* Read in Nsamps samples */
X      if (Scan==EOF || Scan==0)            /*   unless read error or EOF */
X         break;                            /*   in which case we exit and */
X      *Data++ = val;
X      }
X   if (Nsamps > 0)
X      {
X      val = Data - DataStart;              /* (Calc return value) */
X      while (--Nsamps >= 0)                /*   fill unread spaces with 0's */
X         *Data++ = 0;                      /*   and return FALSE */
X      return(val);
X      }
X   return(0);
X}
X
X
X
XwriteData(Nout, Data)
Xint Nout;
XHWORD Data[];
X{
X   short s;
X   /* Write Nout samples to ascii file */
X   while (--Nout >= 0) {
X      s = *Data++;
X      fwrite(&s, 1, 2, fout);
X   }
X}
X
X
X
X
Xgetparams(Factor, InterpFilt, Nmult)
Xdouble *Factor;
XBOOL *InterpFilt;
XUHWORD *Nmult;
X{
X   /* Check for illegal constants */
X   if (Np >= 16)
X      fail("Error: Np>=16");
X   if (Nb+Nhg+NLpScl >= 32)
X      fail("Error: Nb+Nhg+NLpScl>=32");
X   if (Nh+Nb > 32)
X      fail("Error: Nh+Nb>32");
X
X   /* Default conversion parameters */
X   *Factor = 2;
X   *InterpFilt = TRUE;
X
X   /* Set up conversion parameters */
X   while (TRUE)
X      {
X      *Factor =
X         GetDouble("Sampling-rate conversion factor",*Factor,FactorHelp);
X      if (*Factor <= MAXFACTOR)
X         break;
X      printf("Error: Factor must be <= %g to prevent out buffer overflow",
X         MAXFACTOR);
X      *Factor = MAXFACTOR;
X      }
X   *Nmult = GetUHWORD("Nmult",*Nmult,NmultHelp);
X   *InterpFilt = Query("Interpolate filter?", *InterpFilt, InterpFiltHelp);
X}
X
X
X
X/* Sampling rate up-conversion only subroutine;
X * Slightly faster than down-conversion;
X */
XSrcUp(X, Y, Factor, Time, Nx, Nwing, LpScl, Imp, ImpD, Interp)
XHWORD X[], Y[], Imp[], ImpD[];
XUHWORD Nx, Nwing, LpScl;
XUWORD *Time;
Xdouble Factor;
XBOOL Interp;
X{
X   HWORD *Xp, *Ystart;
X   WORD v;
X
X   double dt;                  /* Step through input signal */ 
X   UWORD dtb;                  /* Fixed-point version of Dt */
X   UWORD endTime;              /* When Time reaches EndTime, return to user */
X
X   dt = 1.0/Factor;            /* Output sampling period */
X   dtb = dt*(1<<Np) + 0.5;     /* Fixed-point representation */
X
X   Ystart = Y;
X   endTime = *Time + (1<<Np)*(WORD)Nx;
X   while (*Time < endTime)
X      {
X      Xp = &X[*Time>>Np];      /* Ptr to current input sample */
X      v = FilterUp(Imp, ImpD, Nwing, Interp, Xp, (HWORD)(*Time&Pmask),
X         -1);                  /* Perform left-wing inner product */
X      v += FilterUp(Imp, ImpD, Nwing, Interp, Xp+1, (HWORD)((-*Time)&Pmask),
X         1);                   /* Perform right-wing inner product */
X      v >>= Nhg;               /* Make guard bits */
X      v *= LpScl;              /* Normalize for unity filter gain */
X      *Y++ = v>>NLpScl;        /* Deposit output */
X      *Time += dtb;            /* Move to next sample by time increment */
X      }
X   return (Y - Ystart);        /* Return the number of output samples */
X}
X
X
X
X/* Sampling rate conversion subroutine */
X
XSrcUD(X, Y, Factor, Time, Nx, Nwing, LpScl, Imp, ImpD, Interp)
XHWORD X[], Y[], Imp[], ImpD[];
XUHWORD Nx, Nwing, LpScl;
XUWORD *Time;
Xdouble Factor;
XBOOL Interp;
X{
X   HWORD *Xp, *Ystart;
X   WORD v;
X
X   double dh;                  /* Step through filter impulse response */
X   double dt;                  /* Step through input signal */
X   UWORD endTime;              /* When Time reaches EndTime, return to user */
X   UWORD dhb, dtb;             /* Fixed-point versions of Dh,Dt */
X
X   dt = 1.0/Factor;            /* Output sampling period */
X   dtb = dt*(1<<Np) + 0.5;     /* Fixed-point representation */
X
X   dh = MIN(Npc, Factor*Npc);  /* Filter sampling period */
X   dhb = dh*(1<<Na) + 0.5;     /* Fixed-point representation */
X
X   Ystart = Y;
X   endTime = *Time + (1<<Np)*(WORD)Nx;
X   while (*Time < endTime)
X      {
X      Xp = &X[*Time>>Np];      /* Ptr to current input sample */
X      v = FilterUD(Imp, ImpD, Nwing, Interp, Xp, (HWORD)(*Time&Pmask),
X          -1, dhb);            /* Perform left-wing inner product */
X      v += FilterUD(Imp, ImpD, Nwing, Interp, Xp+1, (HWORD)((-*Time)&Pmask),
X           1, dhb);            /* Perform right-wing inner product */
X      v >>= Nhg;               /* Make guard bits */
X      v *= LpScl;              /* Normalize for unity filter gain */
X      *Y++ = v>>NLpScl;        /* Deposit output */
X      *Time += dtb;            /* Move to next sample by time increment */
X      }
X   return (Y - Ystart);        /* Return the number of output samples */
X}
X
X
X
Xint resample(Imp, ImpD, LpScl, Nmult, Nwing, Factor, InterpFilt)
XHWORD Imp[], ImpD[];
XUHWORD LpScl, Nmult, Nwing;
Xdouble Factor;
XBOOL InterpFilt;
X{
X   UWORD Time;                       /* Current time/pos in input sample */
X   UHWORD Xp, Ncreep, Xoff, Xread;
X   HWORD X[IBUFFSIZE], Y[OBUFFSIZE]; /* I/O buffers */
X   UHWORD Nout, Nx;
X   int i, Ycount, last;
X
X   /* Account for increased filter gain when using Factors less than 1 */
X   if (Factor < 1)
X      LpScl = LpScl*Factor + 0.5;
X   /* Calc reach of LP filter wing & give some creeping room */
X   Xoff = ((Nmult+1)/2.0) * MAX(1.0,1.0/Factor) + 10;
X   if (IBUFFSIZE < 2*Xoff)      /* Check input buffer size */
X      fail("IBUFFSIZE (or Factor) is too small");
X   Nx = IBUFFSIZE - 2*Xoff;     /* # of samples to proccess each iteration */
X
X   last = FALSE;                /* Have not read last input sample yet */
X   Ycount = 0;                  /* Current sample and length of output file */
X   Xp = Xoff;                   /* Current "now"-sample pointer for input */
X   Xread = Xoff;                /* Position in input array to read into */
X   Time = (Xoff<<Np);           /* Current-time pointer for converter */
X
X   for (i=0; i<Xoff; X[i++]=0); /* Need Xoff zeros at begining of sample */
X    
X   do {
X      if (!last)                /* If haven't read last sample yet */
X         {
X         last = readData(X, IBUFFSIZE, (int)Xread); /* Fill input buffer */
X         if (last && (last+Xoff<Nx))  /* If last sample has been read... */
X            Nx = last + Xoff;   /* ...calc last sample affected by filter */
X	 }
X      if (Factor >= 1)          /* Resample stuff in input buffer */
X         Ycount += Nout = SrcUp(X,Y,Factor,&Time,Nx,Nwing,LpScl,Imp,ImpD,
X            InterpFilt);        /* SrcUp() is faster if we can use it */
X      else
X         Ycount += Nout = SrcUD(X,Y,Factor,&Time,Nx,Nwing,LpScl,Imp,ImpD,
X            InterpFilt);
X      Time -= (Nx<<Np);           /* Move converter Nx samples back in time */
X      Xp += Nx;                   /* Advance by number of samples processed */
X      Ncreep = (Time>>Np) - Xoff; /* Calc time accumulation in Time */
X      if (Ncreep)
X         {
X         Time -= (Ncreep<<Np);    /* Remove time accumulation */
X         Xp += Ncreep;            /* and add it to read pointer */
X         }
X      for (i=0; i<IBUFFSIZE-Xp+Xoff; i++)   /* Copy portion of input signal */
X         X[i] = X[i+Xp-Xoff];               /* That must be re-used */
X      if (last)                             /* If near end of sample... */
X         {
X         last -= Xp;             /* ...keep track were it ends */
X         if (!last)              /* Lengthen input by 1 sample if... */
X            last++;              /* ...needed to keep flag TRUE */
X	 }
X      Xread = i;                 /* Pos in input buff to read new data into */
X      Xp = Xoff;
X
X      if (Nout > OBUFFSIZE)      /* Check to see if output buff overflowed */
X         fail("Output array overflow");
X
X      writeData((int)Nout ,Y);   /* Write data in output buff to file */
X      } while (last >= 0);       /* Continue until done processing samples */
X   return(Ycount);               /* Return # of samples in output file */
X}
X
X
X
X
Xmain(argc, argv)
Xint argc;
Xchar *argv[];
X{
X   double Factor;               /* Factor = Fout/Fin */
X   BOOL InterpFilt;             /* TRUE means interpolate filter coeffs */
X   UHWORD LpScl, Nmult, Nwing;  
X   HWORD Imp[MAXNWING];         /* Filter coefficients */
X   HWORD ImpD[MAXNWING];        /* ImpD[n] = Imp[n+1]-Imp[n] */
X   int outCount;
X
X   Nmult = 13;
X   printf("\nSampling Rate Conversion Program\n");
X   printf("Written by:  Julius O. Smith (CCRMA)\n");
X   printf("Translated to C by: Christopher Lee Fraley (cf0v@spice.cs.cmu.edu)\n");
X   printf("Version %s\n", resampleVERSION);
X   openData(argc, argv);         /* Interp arguments and open i&o files */
X   getparams(&Factor, &InterpFilt, &Nmult);
X   if (readFilter(Imp, ImpD, &LpScl, Nmult, &Nwing))
X      fail("could not find filter file, or syntax error in filter file");
X   printf("\nStarting Conversion...\n");
X   outCount = resample(Imp, ImpD, LpScl, Nmult, Nwing, Factor, InterpFilt);
X   closeData();
X
X   printf("...Conversion Finished:  %d output samples.\n\n",outCount);
X}
X
SHAR_EOF
chmod 0644 resample.c ||
echo 'restore of resample.c failed'
Wc_c="`wc -c < 'resample.c'`"
test 13378 -eq "$Wc_c" ||
	echo 'resample.c: original size 13378, current size' "$Wc_c"
fi
# ============= warp.c ==============
if test -f 'warp.c' -a X"$1" != X"-c"; then
	echo 'x - skipping warp.c (File already exists)'
else
echo 'x - extracting warp.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'warp.c' &&
X
X/*
X * FILE: warp.c
X *   BY: Christopher Lee Fraley (cf0v@spice.cs.cmu.edu)
X * NOTE: Based upon SAIL program by Julius O. Smith (CCRMA, Stanford U)
X * DATE: 17-JUN-88
X */
X
Xchar warpVERSION[] = "1.0  24-JUN-88, 4:20pm";
X
X/*
X *      The original SAIL program on which this is based may be found in
X *      /../s/usr/mhs/resample  or  /../x/usr/rbd/src/resample
X *
X *      Implement dynamic sampling-rate changes; uses a second file to
X *      indicate where periods should fall.  This program may be used to add
X *      or remove vibrato and micro pitch variations from a sound sample.
X *      Based on SRCONV.SAI[SYS,MUS] with new algorithm by JOS&PG.
X *      The program internally uses 16-bit data and 16-bit filter
X *      coefficients.
X *
X *      Reference: "A Flexible Sampling-Rate Conversion Method,"
X *      J. O. Smith and P. Gossett, ICASSP, San Diego, 1984, Pgs 19.4.
X *
X *      * Need overflow detection in Filter() routines.  Also, we want
X *        to saturate instead of wrap-around and report number of clipped
X *        samples.
X */
X
X/* CHANGES from original SAIL program:
X *
X * 1. SrcUD() - Switched order of making guard bits (v>>Nhg) and
X *       normalizing.  This was done to prevent overflow.
X */
X
X/* LIBRARIES needed:
X *
X * 1. filterkit
X *       readFilter() - reads standard filter file.
X *       FilterUD()   - applies filter to sample for any Factor.
X *       GetDouble()  - prompt user for a double with help.
X *
X * 2. math
X */
X
X
X
X
X#include <stdio.h>
X#include <math.h>
X#include <string.h>
X#include "stdefs.h"
X#include "filterkit.h"
X#include "resample.h"
X
X#define IBUFFSIZE 1024                         /* Input buffer size */
X#define OBUFFSIZE (IBUFFSIZE*MAXFACTOR+2)      /* Calc'd out buffer size */
X
X
Xfail(s)
Xchar *s;
X{
X   printf("warp: %s\n\n", s);     /* Display error message  */
X   exit(1);                       /* Exit, indicating error */
X}
X
Xfails(s,s2)
Xchar *s, *s2;
X{
X   printf("warp: ");              /* Display error message  */
X   printf(s, s2);
X   printf("\n\n");
X   exit(1);                       /* Exit, indicating error */
X}
X
X
X
X/* Global file pointers: */
XFILE *fin, *fout;
X
X
Xint openData(argc, argv)
Xint argc;
Xchar *argv[];
X{
X   if (argc != 3)
X      fail("format is 'warp <file-in> <file-out>'");
X   if (NULL == (fin = fopen(*++argv,"r")))
X      fails("could not open <file-in> file '%s'",*argv);
X   if (NULL == (fout = fopen(*++argv,"w")))
X      fails("could not open <file-out> file '%s'",*argv);
X}
X
X
XcloseData()
X{
X   (void) fclose(fin);
X   (void) fclose(fout);
X}
X
X
X
Xint readData(Data, DataArraySize, Xoff)  /* return: 0 - notDone */
XHWORD Data[];                            /*        >0 - index of last sample */
Xint DataArraySize, Xoff;
X{
X   int Nsamps, Scan, val;
X   HWORD *DataStart;
X
X   DataStart = Data;
X   Nsamps = DataArraySize - Xoff;   /* Calculate number of samples to get */
X   Data += Xoff;                    /* Start at designated sample number */
X   for (; Nsamps>0; Nsamps--)
X      {
X      Scan = fscanf(fin, "%d", &val);      /* Read in Nsamps samples */
X      if (Scan==EOF || Scan==0)            /*   unless read error or EOF */
X         break;                            /*   in which case we exit and */
X      *Data++ = val;
X      }
X   if (Nsamps > 0)
X      {
X      val = Data - DataStart;              /* (Calc return value) */
X      while (--Nsamps >= 0)                /*   fill unread spaces with 0's */
X         *Data++ = 0;                      /*   and return FALSE */
X      return(val);
X      }
X   return(0);
X}
X
X
X
XwriteData(Nout, Data)
Xint Nout;
XHWORD Data[];
X{
X   while (--Nout >= 0)                  /* Write Nout samples to ascii file */
X      fprintf(fout, "%d\n", *Data++);
X}
X
X
X
X
Xcheck()
X{
X   /* Check for illegal constants */
X   if (Np >= 16)
X      fail("Error: Np>=16");
X   if (Nb+Nhg+NLpScl >= 32)
X      fail("Error: Nb+Nhg+NLpScl>=32");
X   if (Nh+Nb > 32)
X      fail("Error: Nh+Nb>32");
X}
X
X
X/* Globals for warping frequency */
Xdouble a,b,c,d,e,f,Total;
Xint Acc;
X
XinitWarp()
X{
X   Total = GetDouble("\n# of input samples",12000.0,"");
X
X   printf("Warping function:  Fout/Fin = w(n)\n");
X   printf("  w(n) = off + [amp * sin(ph1+frq1*n/N) * sin(ph2+frq2*n/N)]\n");
X   printf("  where: off,amp,ph1   - parameters\n");
X   printf("         frq1,ph2,frq2 - more parameters\n");
X   printf("         n             - current input sample number\n");
X   printf("         N             - total number of input samples\n");
X
X   a = GetDouble("off",1.5,"");
X   b = GetDouble("amp",-0.5,"");
X   c = D2R * GetDouble("ph1 (degrees)",-90.0,"");
X   d = GetDouble("frq1",1.0,"");
X   e = D2R * GetDouble("ph2 (degrees)",90.0,"");
X   f = GetDouble("frq2",0.0,"");
X}
X
X
Xdouble warpFunction(Time)
XUWORD Time;
X{
X   double fTime,t;
X
X   /* Calc absolute position in input file: */
X   fTime = (double)Time / (double)(1<<Np) + (double)Acc;
X   /* Calc fraction of file processed: */
X   t = fTime/Total;
X   /* Return warp factor: */
X   return (1.0 / (a + b*sin(c+PI2*d*t)*sin(e+PI2*f*t)));
X}
X
X
X/* Sampling rate conversion subroutine */
X
XSrcUD(X, Y, Time, Nx, Nwing, LpScl, Imp, ImpD, Interp)
XHWORD X[], Y[], Imp[], ImpD[];
XUHWORD Nx, Nwing, LpScl;
XUWORD *Time;
XBOOL Interp;
X{
X   HWORD *Xp, *Ystart;
X   WORD v;
X
X   UHWORD NewScl;              /* Unity gain scale factor */
X   double dh;                  /* Step through filter impulse response */
X   double dt;                  /* Step through input signal */
X   UWORD endTime;              /* When Time reaches EndTime, return to user */
X   UWORD dhb, dtb;             /* Fixed-point versions of Dh,Dt */
X   double Factor;              /* Current resampling factor */
X
X   Ystart = Y;
X   endTime = *Time + (1<<Np)*(WORD)Nx;
X   while (*Time < endTime)
X      {
X      Factor = warpFunction(*Time);     /* Get new conversion Factor */
X      NewScl = LpScl * Factor;          /* Calc new normalizing scalar */
X      dt = 1.0 / Factor;                /* New output sampling period */
X      dtb = dt*(1<<Np) + 0.5;           /* Fixed-point representation */
X      dh = MIN(Npc, Factor * Npc);      /* New filter sampling period */
X      dhb = dh*(1<<Na) + 0.5;           /* Fixed-point representation */
X      Xp = &X[*Time>>Np];               /* Ptr to current input sample */
X      v = FilterUD(Imp, ImpD, Nwing, Interp, Xp, (HWORD)(*Time&Pmask),
X          -1, dhb);                     /* Perform left-wing inner product */
X      v += FilterUD(Imp, ImpD, Nwing, Interp, Xp+1, (HWORD)((-*Time)&Pmask),
X           1, dhb);                     /* Perform right-wing inner product */
X      v >>= Nhg;                        /* Make guard bits */
X      v *= NewScl;                      /* Normalize for unity filter gain */
X      *Y++ = v>>NLpScl;                 /* Deposit output */
X      *Time += dtb;                     /* Move to next sample by time inc */
X      }
X   return (Y - Ystart);                 /* Return the # of output samples */
X}
X
X
X
Xint resample(Imp, ImpD, LpScl, Nmult, Nwing, InterpFilt)
XHWORD Imp[], ImpD[];
XUHWORD LpScl, Nmult, Nwing;
XBOOL InterpFilt;
X{
X   UWORD Time;                       /* Current time/pos in input sample */
X   UHWORD Xp, Xread, Ncreep, Xoff;
X   HWORD X[IBUFFSIZE], Y[OBUFFSIZE]; /* I/O buffers */
X   UHWORD Nout, Nx;
X   int i, Ycount, last;
X
X   /* Calc reach of LP filter wing & give some creeping room */
X   Xoff = ((Nmult+1)/2.0) * MAX(1.0,1.0*MAXFACTOR) + 10;
X   if (IBUFFSIZE < 2*Xoff)      /* Check input buffer size */
X      fail("IBUFFSIZE (or Factor) is too small");
X   Nx = IBUFFSIZE - 2*Xoff;     /* # of samples to proccess each iteration */
X
X   last = FALSE;                /* Have we read last input sample yet? */
X   Ycount = 0;                  /* Output sample # and length of out file */
X   Xp = Xoff;                   /* Current position in input buffer */
X   Xread = Xoff;                /* Position in input array to read into */
X   Time = (Xoff<<Np);           /* Current-time pointer for converter */
X   Acc = -Xoff;                 /* Acc + int(Time) = index into input file */
X
X   for (i=0; i<Xoff; X[i++]=0); /* Need Xoff zeros at begining of sample */
X
X   do {
X      if (!last)                /* If haven't read last sample yet */
X         {
X         last = readData(X, IBUFFSIZE, (int)Xread); /* Fill input buffer */
X         if (last && (last+Xoff<Nx))  /* If last sample has been read... */
X            Nx = last + Xoff;   /* ...calc last sample affected by filter */
X	 }
X      Ycount += Nout = SrcUD(X,Y,&Time,Nx,Nwing,LpScl,Imp,ImpD,InterpFilt);
X      Time -= (Nx<<Np);           /* Move converter Nx samples back in time */
X      Xp += Nx;                   /* Advance by number of samples processed */
X      Acc += Nx;                  /* We moved Nx samples in the input file */
X      Ncreep = (Time>>Np) - Xoff; /* Calc time accumulation in Time */
X      if (Ncreep)
X         {
X         Time -= (Ncreep<<Np);    /* Remove time accumulation */
X         Xp += Ncreep;            /* and add it to read pointer */
X         Acc += Ncreep;
X         }
X      for (i=0; i<IBUFFSIZE-Xp+Xoff; i++)   /* Copy portion of input signal */
X         X[i] = X[i+Xp-Xoff];               /*    that must be re-used */
X      if (last)                             /* If near end of sample... */
X         {
X         last -= Xp;             /* ...keep track were it ends */
X         if (!last)              /* Lengthen input by 1 sample */
X            last++;              /*  if needed to keep flag TRUE */
X	 }
X      Xread = i;                 /* Pos in input buff to read new data into */
X      Xp = Xoff;
X
X      if (Nout > OBUFFSIZE)      /* Check to see if output buff overflowed */
X         fail("Output array overflow");
X
X      writeData((int)Nout, Y);   /* Write data in output buff to file */
X      } while (last >= 0);       /* Continue until done processing samples */
X   return(Ycount);               /* Return # of samples in output file */
X}
X
X
X
X
Xmain(argc,argv)
Xint argc;
Xchar *argv[];
X{
X   BOOL InterpFilt;             /* TRUE means interpolate filter coeffs */
X   UHWORD LpScl, Nmult, Nwing;
X   HWORD Imp[MAXNWING];         /* Filter coefficients */
X   HWORD ImpD[MAXNWING];        /* ImpD[n] = Imp[n+1]-Imp[n] */
X   int outCount;
X
X   Nmult = 13;
X   InterpFilt = TRUE;
X   printf("\nDynamic Rate Resampling Program (for interesting effects)\n");
X   printf("Written by:  Christopher Lee Fraley (cf0v@spice.cs.cmu.edu)\n");
X   printf("Based upon SAIL program by Julius O. Smith (CCRMA)\n");
X   printf("Version %s\n", warpVERSION);
X   check();                      /* Check constants for validity */
X   openData(argc, argv);         /* Interp arguments and open i&o files */
X   initWarp();                   /* Set up the warp function's parameters */
X   if (readFilter(Imp, ImpD, &LpScl, Nmult, &Nwing))
X      fail("could not open Filter file, or syntax error in filter file");
X   printf("\nWarp Speed Full Ahead...\n");
X   outCount = resample(Imp, ImpD, LpScl, Nmult, Nwing, InterpFilt);
X   closeData();
X
X   printf("...Returning To Ion Drive:  %d output samples.\n\n", outCount);
X}
X
SHAR_EOF
chmod 0644 warp.c ||
echo 'restore of warp.c failed'
Wc_c="`wc -c < 'warp.c'`"
test 10961 -eq "$Wc_c" ||
	echo 'warp.c: original size 10961, current size' "$Wc_c"
fi
exit 0