#!/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
X#include
X#include
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 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= 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>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<>= Nhg;
X v *= LpScl;
X out = (double)v / (double)(1<= 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
X#include
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 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 '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
X#include
X#include
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 '");
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]; /* 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]; /* 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<= 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) - Xoff; /* Calc time accumulation in Time */
X if (Ncreep)
X {
X Time -= (Ncreep< 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
X#include
X#include
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 '");
X if (NULL == (fin = fopen(*++argv,"r")))
X fails("could not open file '%s'",*argv);
X if (NULL == (fout = fopen(*++argv,"w")))
X fails("could not open 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]; /* 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) - Xoff; /* Calc time accumulation in Time */
X if (Ncreep)
X {
X Time -= (Ncreep< 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