Source Code 1994 March/Source_Code_CD-ROM_Walnut_Creek_March_1994.iso/compsrcs/unix/volume26/filtrkit/part01

Text File | 1993-04-14 | 87.6 KB | 2,602 lines

Newsgroups: comp.sources.unix
From: [e-mail] (Christopher Lee Fraley)
Subject: v26i166: filterkit - a collection of audio filters and resamplers, Part01/01
Sender: [e-mail]
Approved: [e-mail]

Submitted-By: [e-mail] (Christopher Lee Fraley)
Posting-Number: Volume 26, Issue 166
Archive-Name: filterkit/part01

     The original SAIL resampling program is now located among four files:
"filterkit.c", "makefilter.c", "resample.c", and "resample.h".  These files
contain:

    filterkit.c - The source code for the library.  Contains routines
        to calculate filter coeff's, convert these to a 16-bit
        integer format, write a filter file, read a filter file,
        apply the filter when up-sampling, apply the filter when
        up-/down-sampling, and apply the filter to find a zero-
        crossing.

    makefilter.c - Routines to design and write a filter to a file.
        Creates a Kaiser-windowed lowpass filter based upon routines
        found in the filterkit library.

    resample.c - The original resampling program stripped down to the
            part that does the actual resampling.  Designing a filter
            was moved to "makefilter.c", the conversion constants
        were moved into "resample.c", and various general filter
        routines were moved into the filterkit library.

    resample.h - Contains the resampling conversion constants.  These
        constants govern such things as the number of bits in the
        input sample & filter coeff's, the number of bits to the
        right of the binary-point for fixed-point math, etc.

     Other programs included in this package are the Makefiles for each of
the programs (and the filterkit library), the needed headerfiles, and a
human-readable copy of the lint version of the filterkit library:

    warp.c - Very similar to "resample.c", except is modified to allow
        the conversion factor to change on a sample-by-sample basis,
        currently using sinusiods.  This program can be used to add
        vibrato or other pitch-shifting effects to a sample.

    stdefs.h - Contains the definitions for words and halfwords, which
        can change from machine to machine.  Also, useful constants
        are defined, as well as a few useful macros.

    filterkit.h - Contains the declarations of each of the functions in
        the filterkit library.  Must be included by any program using
        any of the filterkit routines.

    llib-lfilterkit - A human-readble version of the lint-library for
        filterkit.  Simply defines the type of each function and
        the routine's formal parameters, so lint can check programs
        which use the filterkit library.  The compact version of this
        file is created via the -C option of lint.  See the filterkit
        Makefile for details.

     For further information on any of the files, look in the comments and
help strings contained within that file.  Also I can be contacted at:
    [e-mail], or
    [e-mail].

     --Christopher Lee Fraley

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 1 (of 1)."
# Contents:  MANIFEST Makefile Notes README.filterkit README.resample
#   filterkit.c filterkit.h kaiser.c makefilter.c protos.h resample.c
#   resample.h stdefs.h warp.c
# Wrapped by [e-mail] on Thu Apr 15 18:45:31 1993
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'MANIFEST' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'MANIFEST'\"
else
echo shar: Extracting \"'MANIFEST'\" \(548 characters\)
sed "s/^X//" >'MANIFEST' <<'END_OF_FILE'
X   File Name        Archive #    Description
X-----------------------------------------------------------
X MANIFEST                   1    This shipping list
X Makefile                   1
X Notes                      1
X README.filterkit           1
X README.resample            1
X filterkit.c                1
X filterkit.h                1
X kaiser.c                   1
X makefilter.c               1
X protos.h                   1
X resample.c                 1
X resample.h                 1
X stdefs.h                   1
X warp.c                     1
END_OF_FILE
if test 548 -ne `wc -c <'MANIFEST'`; then
    echo shar: \"'MANIFEST'\" unpacked with wrong size!
fi
# end of 'MANIFEST'
fi
if test -f 'Makefile' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Makefile'\"
else
echo shar: Extracting \"'Makefile'\" \(954 characters\)
sed "s/^X//" >'Makefile' <<'END_OF_FILE'
SRC = kaiser.c filterkit.c makefilter.c resample.c warp.c
OBJ = kaiser.o filterkit.o makefilter.o resample.o warp.o
BIN = kaiser makefilter resample warp
X
XFILES = Readme Notes Makefile filterkit.h protos.h resample.h stdefs.h $(SRC)
X
INSTALLDIR = /usr/local/bin
X
all: $(BIN)
X
install: $(BIN)
X    cp $(BIN) "$(INSTALLDIR)"
X
kaiser: kaiser.c filterkit.c
X    rm -f kaiser
X    $(CC) $(CFLAGS) kaiser.c filterkit.c -o kaiser -lm
X
makefilter: makefilter.c filterkit.c
X    rm -f makefilter
X    $(CC) $(CFLAGS) makefilter.c filterkit.c -o makefilter -lm
X
resample: resample.c filterkit.c
X    rm -f resample
X    $(CC) $(CFLAGS) resample.c filterkit.c -o resample -lm
X
warp: warp.c filterkit.c
X    rm -f warp
X    $(CC) $(CFLAGS) warp.c filterkit.c -o warp -lm
X
X$(SRC): filterkit.h stdefs.h
resample.c: resample.h
X
X# Shar: -F (prefix all lines with X),
X#     -s addr (set return addr of poster)
shar: $(FILES)
X    /usr2/tools/shar/shar -F -l 50 -o shar -n resample -s [e-mail] $(FILES)
END_OF_FILE
if test 954 -ne `wc -c <'Makefile'`; then
    echo shar: \"'Makefile'\" unpacked with wrong size!
fi
# end of 'Makefile'
fi
if test -f 'Notes' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Notes'\"
else
echo shar: Extracting \"'Notes'\" \(9614 characters\)
sed "s/^X//" >'Notes' <<'END_OF_FILE'
X
X                NOTES
X
X
PROGRAM ORGANIZATION:
X
X     The original SAIL resampling program is now located among four files:
X"filterkit.c", "makefilter.c", "resample.c", and "resample.h".  These files
contain:
X
X    filterkit.c - The source code for the library.  Contains routines
X        to calculate filter coeff's, convert these to a 16-bit
X        integer format, write a filter file, read a filter file,
X        apply the filter when up-sampling, apply the filter when
X        up-/down-sampling, and apply the filter to find a zero-
X        crossing.
X
X    makefilter.c - Routines to design and write a filter to a file.
X        Creates a Kaiser-windowed lowpass filter based upon routines
X        found in the filterkit library.
X
X    resample.c - The original resampling program stripped down to the
X            part that does the actual resampling.  Designing a filter
X            was moved to "makefilter.c", the conversion constants
X        were moved into "resample.c", and various general filter
X        routines were moved into the filterkit library.
X
X    resample.h - Contains the resampling conversion constants.  These
X        constants govern such things as the number of bits in the
X        input sample & filter coeff's, the number of bits to the
X        right of the binary-point for fixed-point math, etc.
X
X     Other programs included in this package are the Makefiles for each of
the programs (and the filterkit library), the needed headerfiles, and a
human-readable copy of the lint version of the filterkit library:
X
X    warp.c - Very similar to "resample.c", except is modified to allow
X        the conversion factor to change on a sample-by-sample basis,
X        currently using sinusiods.  This program can be used to add
X        vibrato or other pitch-shifting effects to a sample.
X
X    stdefs.h - Contains the definitions for words and halfwords, which
X        can change from machine to machine.  Also, useful constants
X        are defined, as well as a few useful macros.
X
X    filterkit.h - Contains the declarations of each of the functions in
X        the filterkit library.  Must be included by any program using
X        any of the filterkit routines.
X
X    llib-lfilterkit - A human-readble version of the lint-library for
X        filterkit.  Simply defines the type of each function and
X        the routine's formal parameters, so lint can check programs
X        which use the filterkit library.  The compact version of this
X        file is created via the -C option of lint.  See the filterkit
X        Makefile for details.
X
X
X
NOTES ON COMPILATION:
X
X     All of the modules should compile without any problems, although search
paths will have to be set up correctly so "cc" can find all of the proper
include and library files.  Alternatively, one could also explicly cite an
absolute file path for each of the files.
X
X     A few minor changes have been made to some of these files since they were
last compiled:
X    1. All absolute paths were removed from the Makefiles to make
X        adaption to your site less confusing.
X    2. The source for the filterkit library was concatinated into one
X        file.  The Query(), GetUShort(), GetDouble(), and GetString()
X        routines were originally in a separate file "query.c".
X    3. "llib-lfilterkit" was modified to reflect the latest change to
X        the zerox() routine.  This has not been tested yet.
X
X     The routines originally in "query.c" may have to be modified slightly
to get them to compile.  Apparently, the getstr() routine (which is used in
each of these routines) is not a part of the standard C library on all Unix
systems.  Its specification is:
X    getstr(prompt, defalt, answer)
X    char *prompt, *defalt, *answer;
X    Getstr() will print the passed "prompt" as a message to the user, and
X        wait for the user to type an input string.  The string is
X        then copied into "answer".  If the user types just a carriage
X        return, then the string "defalt" is copied into "answer".
X        "Answer" and "defalt" may be the same string, in which case,
X        the default value will be the original contents of "answer".
X
X     Also note that the constants in "resample.h" and the typedefs in
X"stdefs.h" may need to be changed to reflect the word-size of the machine
the program will be running on.  These files currently reflect a word-size
of 32-bits.
X
X
X
STRUCTURE OF FILES:
X
X     There are two types of files used to pass information from one program
the the next: the filter file and a sound-sample file.  The filter file
has the following format:
X
X    File Name: "F" Nmult "T" Nhc ".filter"
X        example: "F13T8.filter" and "F27T8.filter"
X
X    Structure of File:
X        "ScaleFactor" LpScl
X        "Length" Nwing
X        "Coeffs:"
X        Imp[0]
X        Imp[1]
X          :
X        Imp[Nwing-1]
X        "Differences:"
X        ImpD[0]
X        ImpD[1]
X          :
X        ImpD[Nwing-1]
X        EOF
X
X    where:    Something enclosed in "" inicates specific characters found
X            in the file.
X        Nmult, Nwing, Imp[], and ImpD[] are variables (HWORD)
X        Npc is a conversion constant.
X        EOF is the end of the file.
X    See writeFilter() and readFilter() in "filterkit.c" for more details.
X
X
The sound-sample file has the following format:
X
X    File Name: unspecified.
X        example: "horn.in", "Cello.adc", etc.
X
X    Structure of File:
X        Data[0]
X        Data[1]
X          :
X        Data[n]
X        EOF
X
X    where:    Data[i] is a 16-bit integer in ascii format.  (Example:
X            "32767", "-232")
X        EOF is the end of the file.
X        Note that there is no indication of the length of the
X            sample;  the EOF is the only indication that the
X            last sample has been reached.
X
X
X
SUGGESTIONS, CORRECTIONS, IMPROVEMENTS:
X
X     The following were assumed to be mistakes in the original SAIL program,
and have been changed in "resample.c" and other related programs:
X
X    1. LpScl was scaled by Factor (when Factor<1.0) only when the filter
X        was loaded from a file, but not when the filter was created
X        at run-time.
X    2. ImpD[Nwing-1] was originally set to Imp[Nwing-1], but in general
X        ImpD[i] = ImpD[i+1] - ImpD[i], indicating ImpD[Nwing-1]
X        should be set to NEGATIVE Imp[Nwing-1].
X    3. To keep the # of multiplies consistant (and odd), I added a check
X        near the begining of the Filter() routines to skip the last
X        coeff in the right wing when the phase was 180 degrees.
X        (The "if (Inc == 1) End--;" clause).
X    4. When the phase was zero, and the right-wing Filter() is performed,
X        the pointers to the filter coeffs must be incremented to
X        keep the pointing to the correct coeff.  (The "if (Ph == 0)
X        { Hp += Npc; Hdp += Npc; }" clause).
X
X
X     The following were changes made to the basic implementation of the
algorithm:
X
X    1. ImpD[] is now created from the integer Imp[] instead of the
X        double ImpR[], in order to avoid problems with rounding.
X    2. In SrcUp() and SrcUD(), the order of making guard bits (v>>Nhg)
X        and normalizing was switched.  This was done to prevent
X        overflow problems.
X    3. Rounding was removed from SrcUp() and SrcUD() to get more accurate
X        response to a DC input.
X
X
X     The following are improvements that I recommend be made to various
routines:
X
X    1. Change file format of sound files (and possibly filter files) to
X        a binary format.  If you are in a Unix environment, you can
X        use the functions read() and write() to perform low-level
X        i/o, instead of fprintf() and fscanf().  This would be MUCH
X        faster, but I have not made the changes yet because the
X        current code is much clearer, and probably easier to port
X        to your system and your sound-file format.
X    2. I believe the restriction that Nmult must be odd may be removed.
X    3. When Filter() is currently called for the right-wing, Xp+1 must
X        be passed to it.  This should probably be changed so that
X        Filter() adds one to Xp if it is doing the right-wing,
X        instead of the user.  This way, that increment would be
X        transparent to the user, and provide a more consistent
X        interface to these functions.
X
X
X
NOTES ON SPECIFIC FILES:
X
X     "resample.c" - should be a faithful reproduction of the original SAIL
X    program in C.
X
X     "warp.c" - I modified SrcUD() to get the current value for Factor from
X    the function warpFunction().  WarpFunction() should be changed to
X    some form desirable.
X
X     "filterkit.c" - Contains the following useful routines:
X    LpFilter() - Calculates the filter coeffs for a Kaiser-windowed
X        low-pass filter with a given roll-off frequency.  These
X        coeffs are stored into a array of doubles.
X    writeFilter() - Writes a filter to a file.
X    makeFilter() - A section of the original SAIL program.  Calls
X        LpFilter() to create a filter, then scales the double
X        coeffs into a 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.  Both are repoductions of the original SAIL
X        program.
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
XFURTHER INFORMATION:
X
X     Any printf()'s that start with a "|" as the first character in the
control string are used for debugging to to alert the user of an error
condition.
X     For further information on any of the files, look in the comments and
help strings contained within that file.  Also I can be contacted at:
X    [e-mail], or
X    [e-mail].
X
X     --Christopher Lee Fraley
X
END_OF_FILE
if test 9614 -ne `wc -c <'Notes'`; then
    echo shar: \"'Notes'\" unpacked with wrong size!
fi
# end of 'Notes'
fi
if test -f 'README.filterkit' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'README.filterkit'\"
else
echo shar: Extracting \"'README.filterkit'\" \(459 characters\)
sed "s/^X//" >'README.filterkit' <<'END_OF_FILE'
X
X
XFilterkit is a C translation of an old SAIL program (SAIL =
Stanford AI Lab) which does FIR filtering with sample rate
changing on the fly.
X
Christopher Fraley @CMU did the translation and split it up
into a toolkit.  Being ignorant, I've glommed some of the pieces
together into a program that performs straight filtering on
a sound sample file of signed words.
X
The warp program in particular looks very promising; I haven't
used it yet.
X
X    Lance Norskog
END_OF_FILE
if test 459 -ne `wc -c <'README.filterkit'`; then
    echo shar: \"'README.filterkit'\" unpacked with wrong size!
fi
# end of 'README.filterkit'
fi
if test -f 'README.resample' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'README.resample'\"
else
echo shar: Extracting \"'README.resample'\" \(1605 characters\)
sed "s/^X//" >'README.resample' <<'END_OF_FILE'
X(Message inbox:3096)
Received: from netcom4.netcom.com by charon.cwi.nl with SMTP
X    id AA27402 (5.65b/3.8/CWI-Amsterdam); Mon, 8 Mar 1993 22:46:14 +0100
Received: by netcom4.netcom.com (5.65/SMI-4.1/Netcom)
X    id AA04106; Mon, 8 Mar 93 13:45:12 -0800
Date: Mon, 8 Mar 93 13:45:12 -0800
XFrom: [e-mail] (Technically Sweet)
Message-Id: <[e-mail]>
In-Reply-To: [e-mail]
X       "Re: SOX wav.c patch" (Mar  8, 10:07pm)
XX-Mailer: Mail User's Shell (7.2.3 5/22/91)
To: [e-mail]
Subject: resampler part 1 of 2
X
Here it is.  Now I realise there's no help.
X
X'kaiser' reads and writes signed 16-bit samples, but
doesn't care about sampling rates.
X
X'kaiser 2.0' doubles the sampling rate.
X'kaiser 0.5 [ 0.5 ]' cuts the sampling rate in 1/2.
X'kaiser 0.5 0.7' does the same, but with a more gradual slope.
X
I checked this out with the FFT analyzer in MixVIew,
and it really does better than interpolating.
Interp upward reflects the base sound in the higher
frequencies.  This doesn't.
X
X(Message inbox:3097)
Received: from netcom4.netcom.com by charon.cwi.nl with SMTP
X    id AA27418 (5.65b/3.8/CWI-Amsterdam); Mon, 8 Mar 1993 22:46:30 +0100
Received: by netcom4.netcom.com (5.65/SMI-4.1/Netcom)
X    id AA04154; Mon, 8 Mar 93 13:45:29 -0800
Date: Mon, 8 Mar 93 13:45:29 -0800
XFrom: [e-mail] (Technically Sweet)
Message-Id: <[e-mail]>
In-Reply-To: [e-mail]
X       "Re: SOX wav.c patch" (Mar  8, 10:07pm)
XX-Mailer: Mail User's Shell (7.2.3 5/22/91)
To: [e-mail]
Subject: resample part 2 of 2
X
END_OF_FILE
if test 1605 -ne `wc -c <'README.resample'`; then
    echo shar: \"'README.resample'\" unpacked with wrong size!
fi
# end of 'README.resample'
fi
if test -f 'filterkit.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'filterkit.c'\"
else
echo shar: Extracting \"'filterkit.c'\" \(16885 characters\)
sed "s/^X//" >'filterkit.c' <<'END_OF_FILE'
X
X/*
X * FILE: filterkit.c (source for library "filterkit.a")
X *   BY: Christopher Lee Fraley ([e-mail])
X * DESC: Makes a Kaiser-windowed low-pass filter.
X * DATE: 7-JUN-88
X */
X
char 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
double Izero(x)
double 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
void LpFilter(c,N,frq,Beta,Num)
double c[], frq, Beta;
int 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
int writeFilter(Imp, ImpD, LpScl, Nmult, Nwing)
HWORD Imp[], ImpD[];
UHWORD 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: */
static double ImpR[MAXNWING];
X
int makeFilter(Imp, ImpD, LpScl, Nwing, Froll, Beta)
HWORD Imp[], ImpD[];
UHWORD *LpScl, Nwing;
double 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 */
int readFilter(Imp, ImpD, LpScl, Nmult, Nwing)
HWORD Imp[], ImpD[];
UHWORD *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
WORD FilterUp(Imp, ImpD, Nwing, Interp, Xp, Ph, Inc)
HWORD Imp[], ImpD[], *Xp, Inc, Ph;
UHWORD Nwing;
BOOL 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
WORD FilterUD(Imp, ImpD, Nwing, Interp, Xp, Ph, Inc, dhb)
HWORD Imp[], ImpD[], *Xp, Ph, Inc;
UHWORD Nwing, dhb;
BOOL 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
static UHWORD LpScl, Nmult, Nwing;
static HWORD Imp[MAXNWING];
static 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 */
int initZerox(tempNmult)
UHWORD 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
double zerox(Data, Factor)
HWORD *Data;
double 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
BOOL Query(prompt, deflt, help)
char *help, *prompt;
BOOL 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
char *GetString(prompt, deflt, help)
char *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
getstr(s1, deflt, s2)
char *s1, *deflt, *s2;
X{
X    write(1, s1, strlen(s1));
X    read(0, s2, 80);
X}
X
X
X
double GetDouble(title, deflt, help)
char *help, *title;
double 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
unsigned short GetUShort(title, deflt, help)
char *help, *title;
unsigned 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
END_OF_FILE
if test 16885 -ne `wc -c <'filterkit.c'`; then
    echo shar: \"'filterkit.c'\" unpacked with wrong size!
fi
# end of 'filterkit.c'
fi
if test -f 'filterkit.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'filterkit.h'\"
else
echo shar: Extracting \"'filterkit.h'\" \(380 characters\)
sed "s/^X//" >'filterkit.h' <<'END_OF_FILE'
X
X/*
X * FILE:filterkit.h
X *   BY: Christopher Lee Fraley ([e-mail])
X * DATE: 17-JUN-88
X * VERS: 1.0  (20-JUN-88, 3:10pm)
X */
X
void LpFilter();
int writeFilter(),  readFilter(), makeFilter();
WORD FilterUp(), FilterUD();
double GetDouble();
unsigned short GetUShort();
char *GetString();
BOOL Query();
double cubic(), zerox();
X
X#define GetUHWORD(x,y,z) GetUShort(x,y,z)
X
END_OF_FILE
if test 380 -ne `wc -c <'filterkit.h'`; then
    echo shar: \"'filterkit.h'\" unpacked with wrong size!
fi
# end of 'filterkit.h'
fi
if test -f 'kaiser.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'kaiser.c'\"
else
echo shar: Extracting \"'kaiser.c'\" \(12170 characters\)
sed "s/^X//" >'kaiser.c' <<'END_OF_FILE'
X
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 *          ([e-mail] or @andrew.cmu.edu)
X * DATE: 7-JUN-88
X *
X * HACKED: Lance Norskog, Dec. 28, 1991, to make sealed resampler for SoundKit
X */
X
char 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#include "protos.h"
X
X#define IBUFFSIZE 1024                         /* Input buffer size */
X#define OBUFFSIZE (IBUFFSIZE*MAXFACTOR+2)      /* Calc'd out buffer size */
X
fail(s)
char *s;
X{
X   fprintf(stderr, "kaiser: %s\n\n", s);  /* Display error message  */
X   exit(1);                        /* Exit, indicating error */
X}
X
fails(s,s2)
char *s, *s2;
X{
X   fprintf(stderr, "kaiser: ");           /* Display error message  */
X   fprintf(stderr, s,s2);
X   fprintf(stderr, "\n\n");
X   exit(1);                        /* Exit, indicating error */
X}
X
X
X
X
X/* Help Declarations */
X
X/* Global file pointers: */
XFILE *fin, *fout;
X
closeData()
X{
X   (void) fclose(fin);
X   (void) fclose(fout);
X}
X
X
int readData(Data, DataArraySize, Xoff)  /* return: 0 - notDone */
HWORD Data[];                            /*        >0 - index of last sample */
int 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
writeData(Nout, Data)
int Nout;
HWORD 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
getparams(argc, argv, Factor, Froll)
int argc;
char **argv;
double *Factor, *Froll;
X{
X   if ((argc != 2) && (argc != 3))
X      fail("format is 'resample factor [ rolloff ]'");
X   if ((*Factor = atof(argv[1])) < 0.01)
X      fail("conversion factor no good");
X   if (argc == 2)
X      *Froll = 0.5;
X   else if ((*Froll = atof(argv[2])) < 0.01)
X      fail("rolloff factor no good");
X   fin = stdin;
X   fout = stdout;
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
X/* Sampling rate up-conversion only subroutine;
X * Slightly faster than down-conversion;
X */
SrcUp(X, Y, Factor, Time, Nx, Nwing, LpScl, Imp, ImpD, Interp)
HWORD X[], Y[], Imp[], ImpD[];
UHWORD Nx, Nwing, LpScl;
UWORD *Time;
double Factor;
BOOL 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
SrcUD(X, Y, Factor, Time, Nx, Nwing, LpScl, Imp, ImpD, Interp)
HWORD X[], Y[], Imp[], ImpD[];
UHWORD Nx, Nwing, LpScl;
UWORD *Time;
double Factor;
BOOL 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
int resample(Imp, ImpD, LpScl, Nmult, Nwing, Factor, InterpFilt)
HWORD Imp[], ImpD[];
UHWORD LpScl, Nmult, Nwing;
double Factor;
BOOL 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
main(argc, argv)
int argc;
char *argv[];
X{
X   double Factor;               /* Factor = Fout/Fin */
X   double Froll;        /* roll-off frequency */
X   BOOL InterpFilt = TRUE;      /* 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 = 37;
X   getparams(argc, argv, &Factor, &Froll);
X   genFilter(Imp, ImpD, &LpScl, Nmult, &Nwing, Froll);
X   resample(Imp, ImpD, LpScl, Nmult, Nwing, Factor, InterpFilt);
X   closeData();
X}
X
char *cantmake[5] = {
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
genFilter(Imp, ImpD, LpScl, Nmult, Nwing, Froll)
HWORD Imp[MAXNWING];               /* Filter coefficients */
HWORD ImpD[MAXNWING];              /* ImpD[i] = ImpD[i+1] - ImpD[i] */
UHWORD Nmult, *LpScl, *Nwing;
double Froll;
X{
X   double Beta = 2.120;
X   int err;
X
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   if ((Froll<=0) || (Froll>1))
X         fail("Error: Roll-off freq must be 0<Factor<=1\n");
X   if (err = makeFilter(Imp, ImpD, LpScl, *Nwing, Froll, Beta))
X         fails("Error: Unable to make filter: %s\n", cantmake[err]);
X}
X
END_OF_FILE
if test 12170 -ne `wc -c <'kaiser.c'`; then
    echo shar: \"'kaiser.c'\" unpacked with wrong size!
fi
# end of 'kaiser.c'
fi
if test -f 'makefilter.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'makefilter.c'\"
else
echo shar: Extracting \"'makefilter.c'\" \(4633 characters\)
sed "s/^X//" >'makefilter.c' <<'END_OF_FILE'
X
X/*
X * FILE: makefilter.c
X *   BY: Christopher Lee Fraley ([e-mail] or [e-mail])
X * DESC: Makes a Kaiser-windowed low-pass filter.
X * DATE: 7-JUN-88
X */
X
char 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
char 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
char 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
char 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
main()
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
END_OF_FILE
if test 4633 -ne `wc -c <'makefilter.c'`; then
    echo shar: \"'makefilter.c'\" unpacked with wrong size!
fi
# end of 'makefilter.c'
fi
if test -f 'protos.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'protos.h'\"
else
echo shar: Extracting \"'protos.h'\" \(3396 characters\)
sed "s/^X//" >'protos.h' <<'END_OF_FILE'
X#if defined(__STDC__) || defined(__cplusplus)
X# define P_(s) s
X#else
X# define P_(s) ()
X#endif
X
X
X/* filterkit.c */
double Izero P_((double x));
void LpFilter P_((double c[], int N, double frq, double Beta, int Num));
int writeFilter P_((HWORD Imp[], HWORD ImpD[], UHWORD LpScl, UHWORD Nmult, UHWORD Nwing));
int makeFilter P_((HWORD Imp[], HWORD ImpD[], UHWORD *LpScl, UHWORD Nwing, double Froll, double Beta));
int readFilter P_((HWORD Imp[], HWORD ImpD[], UHWORD *LpScl, UHWORD Nmult, UHWORD *Nwing));
WORD FilterUp P_((HWORD Imp[], HWORD ImpD[], UHWORD Nwing, BOOL Interp, HWORD *Xp, HWORD Ph, HWORD Inc));
WORD FilterUD P_((HWORD Imp[], HWORD ImpD[], UHWORD Nwing, BOOL Interp, HWORD *Xp, HWORD Ph, HWORD Inc, UHWORD dhb));
int initZerox P_((UHWORD tempNmult));
double zerox P_((HWORD *Data, double Factor));
BOOL Query P_((char *prompt, BOOL deflt, char *help));
char *GetString P_((char *prompt, char *deflt, char *help));
int getstr P_((char *s1, char *deflt, char *s2));
double GetDouble P_((char *title, double deflt, char *help));
unsigned short GetUShort P_((char *title, unsigned int deflt, char *help));
X
X/* makefilter.c */
int main P_((void));
X
X/* nres.c */
int fail P_((char *s));
int fails P_((char *s, char *s2));
int closeData P_((void));
int readData P_((HWORD Data[], int DataArraySize, int Xoff));
int writeData P_((int Nout, HWORD Data[]));
int getparams P_((int argc, char **argv, double *Factor, double *Froll));
int SrcUp P_((HWORD X[], HWORD Y[], double Factor, UWORD *Time, UHWORD Nx, UHWORD Nwing, UHWORD LpScl, HWORD Imp[], HWORD ImpD[], BOOL Interp));
int SrcUD P_((HWORD X[], HWORD Y[], double Factor, UWORD *Time, UHWORD Nx, UHWORD Nwing, UHWORD LpScl, HWORD Imp[], HWORD ImpD[], BOOL Interp));
int resample P_((HWORD Imp[], HWORD ImpD[], UHWORD LpScl, UHWORD Nmult, UHWORD Nwing, double Factor, BOOL InterpFilt));
int main P_((int argc, char *argv[]));
int genFilter P_((HWORD Imp[MAXNWING ], HWORD ImpD[MAXNWING ], UHWORD *LpScl, UHWORD Nmult, UHWORD *Nwing, double Froll));
X
X/* resample.c */
int fail P_((char *s));
int fails P_((char *s, char *s2));
int openData P_((int argc, char *argv[]));
int closeData P_((void));
int readData P_((HWORD Data[], int DataArraySize, int Xoff));
int writeData P_((int Nout, HWORD Data[]));
int getparams P_((double *Factor, BOOL *InterpFilt, UHWORD *Nmult));
int SrcUp P_((HWORD X[], HWORD Y[], double Factor, UWORD *Time, UHWORD Nx, UHWORD Nwing, UHWORD LpScl, HWORD Imp[], HWORD ImpD[], BOOL Interp));
int SrcUD P_((HWORD X[], HWORD Y[], double Factor, UWORD *Time, UHWORD Nx, UHWORD Nwing, UHWORD LpScl, HWORD Imp[], HWORD ImpD[], BOOL Interp));
int resample P_((HWORD Imp[], HWORD ImpD[], UHWORD LpScl, UHWORD Nmult, UHWORD Nwing, double Factor, BOOL InterpFilt));
int main P_((int argc, char *argv[]));
X
X/* warp.c */
int fail P_((char *s));
int fails P_((char *s, char *s2));
int openData P_((int argc, char *argv[]));
int closeData P_((void));
int readData P_((HWORD Data[], int DataArraySize, int Xoff));
int writeData P_((int Nout, HWORD Data[]));
int check P_((void));
int initWarp P_((void));
double warpFunction P_((UWORD Time));
int SrcUD P_((HWORD X[], HWORD Y[], UWORD *Time, UHWORD Nx, UHWORD Nwing, UHWORD LpScl, HWORD Imp[], HWORD ImpD[], BOOL Interp));
int resample P_((HWORD Imp[], HWORD ImpD[], UHWORD LpScl, UHWORD Nmult, UHWORD Nwing, BOOL InterpFilt));
int main P_((int argc, char *argv[]));
X
X#undef P_
END_OF_FILE
if test 3396 -ne `wc -c <'protos.h'`; then
    echo shar: \"'protos.h'\" unpacked with wrong size!
fi
# end of 'protos.h'
fi
if test -f 'resample.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'resample.c'\"
else
echo shar: Extracting \"'resample.c'\" \(13378 characters\)
sed "s/^X//" >'resample.c' <<'END_OF_FILE'
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 *          ([e-mail] or @andrew.cmu.edu)
X * DATE: 7-JUN-88
X */
X
char 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
fail(s)
char *s;
X{
X   printf("resample: %s\n\n", s);  /* Display error message  */
X   exit(1);                        /* Exit, indicating error */
X}
X
fails(s,s2)
char *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 */
char 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
char 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
char 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
openData(argc,argv)
int argc;
char *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
closeData()
X{
X   (void) fclose(fin);
X   (void) fclose(fout);
X}
X
X
int readData(Data, DataArraySize, Xoff)  /* return: 0 - notDone */
HWORD Data[];                            /*        >0 - index of last sample */
int 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
writeData(Nout, Data)
int Nout;
HWORD 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
getparams(Factor, InterpFilt, Nmult)
double *Factor;
BOOL *InterpFilt;
UHWORD *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 */
SrcUp(X, Y, Factor, Time, Nx, Nwing, LpScl, Imp, ImpD, Interp)
HWORD X[], Y[], Imp[], ImpD[];
UHWORD Nx, Nwing, LpScl;
UWORD *Time;
double Factor;
BOOL 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
SrcUD(X, Y, Factor, Time, Nx, Nwing, LpScl, Imp, ImpD, Interp)
HWORD X[], Y[], Imp[], ImpD[];
UHWORD Nx, Nwing, LpScl;
UWORD *Time;
double Factor;
BOOL 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
int resample(Imp, ImpD, LpScl, Nmult, Nwing, Factor, InterpFilt)
HWORD Imp[], ImpD[];
UHWORD LpScl, Nmult, Nwing;
double Factor;
BOOL 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
main(argc, argv)
int argc;
char *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 ([e-mail])\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
END_OF_FILE
if test 13378 -ne `wc -c <'resample.c'`; then
    echo shar: \"'resample.c'\" unpacked with wrong size!
fi
# end of 'resample.c'
fi
if test -f 'resample.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'resample.h'\"
else
echo shar: Extracting \"'resample.h'\" \(2942 characters\)
sed "s/^X//" >'resample.h' <<'END_OF_FILE'
X
X/*
X * FILE: resample.h
X *   BY: Julius Smith (at CCRMA, Stanford U)
X * C BY: translated from SAIL to C by Christopher Lee Fraley
X *          ([e-mail])
X * DATE: 7-JUN-88
X * VERS: 2.0  (17-JUN-88, 3:00pm)
X */
X
X#define MAXNWING  5122
X#define MAXFACTOR 4       /* Maximum Factor without output buff overflow */
X
X
X
X/* Conversion constants */
X#define Nhc       8
X#define Na        7
X#define Np       (Nhc+Na)
X#define Npc      (1<<Nhc)
X#define Amask    ((1<<Na)-1)
X#define Pmask    ((1<<Np)-1)
X#define Nh       16
X#define Nb       16
X#define Nhxn     14
X#define Nhg      (Nh-Nhxn)
X#define NLpScl   13
X
X/* Description of constants:
X *
X * Npc - is the number of look-up values available for the lowpass filter
X *    between the beginning of its impulse response and the "cutoff time"
X *    of the filter.  The cutoff time is defined as the reciprocal of the
X *    lowpass-filter cut off frequence in Hz.  For example, if the
X *    lowpass filter were a sinc function, Npc would be the index of the
X *    impulse-response lookup-table corresponding to the first zero-
X *    crossing of the sinc function.  (The inverse first zero-crossing
X *    time of a sinc function equals its nominal cutoff frequency in Hz.)
X *    Npc must be a power of 2 due to the details of the current
X *    implementation. The default value of 512 is sufficiently high that
X *    using linear interpolation to fill in between the table entries
X *    gives approximately 16-bit accuracy in filter coefficients.
X *
X * Nhc - is log base 2 of Npc.
X *
X * Na - is the number of bits devoted to linear interpolation of the
X *    filter coefficients.
X *
X * Np - is Na + Nhc, the number of bits to the right of the binary point
X *    in the integer "time" variable. To the left of the point, it indexes
X *    the input array (X), and to the right, it is interpreted as a number
X *    between 0 and 1 sample of the input X.  Np must be less than 16 in
X *    this implementation.
X *
X * Nh - is the number of bits in the filter coefficients. The sum of Nh and
X *    the number of bits in the input data (typically 16) cannot exceed 32.
X *    Thus Nh should be 16.  The largest filter coefficient should nearly
X *    fill 16 bits (32767).
X *
X * Nb - is the number of bits in the input data. The sum of Nb and Nh cannot
X *    exceed 32.
X *
X * Nhxn - is the number of bits to right shift after multiplying each input
X *    sample times a filter coefficient. It can be as great as Nh and as
X *    small as 0. Nhxn = Nh-2 gives 2 guard bits in the multiply-add
X *    accumulation.  If Nhxn=0, the accumulation will soon overflow 32 bits.
X *
X * Nhg - is the number of guard bits in mpy-add accumulation (equal to Nh-Nhxn).
X *
X * NLpScl - is the number of bits allocated to the unity-gain normalization
X *    factor.  The output of the lowpass filter is multiplied by LpScl and
X *    then right-shifted NLpScl bits. To avoid overflow, we must have
X *    Nb+Nhg+NLpScl < 32.
X */
X
END_OF_FILE
if test 2942 -ne `wc -c <'resample.h'`; then
    echo shar: \"'resample.h'\" unpacked with wrong size!
fi
# end of 'resample.h'
fi
if test -f 'stdefs.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'stdefs.h'\"
else
echo shar: Extracting \"'stdefs.h'\" \(699 characters\)
sed "s/^X//" >'stdefs.h' <<'END_OF_FILE'
X
X/*
X * FILE: stdefs.h
X *   BY: Christopher Lee Fraley
X * DESC: Defines standard stuff for inclusion in C programs.
X * DATE: 6-JUN-88
X * VERS: 1.0  (6-JUN-88, 2:00pm)
X */
X
X
X#define TRUE  1
X#define FALSE 0
X
X#define PI (3.14159265358979232846)
X#define PI2 (6.28318530717958465692)
X#define D2R (0.01745329348)          /* (2*pi)/360 */
X#define R2D (57.29577951)            /* 360/(2*pi) */
X
X#define MAX(x,y) ((x)>(y) ?(x):(y))
X#define MIN(x,y) ((x)<(y) ?(x):(y))
X#define ABS(x)   ((x)<0   ?(-(x)):(x))
X#define SGN(x)   ((x)<0   ?(-1):((x)==0?(0):(1)))
X
typedef char           BOOL;
typedef short          HWORD;
typedef unsigned short UHWORD;
typedef int            WORD;
typedef unsigned int   UWORD;
X
END_OF_FILE
if test 699 -ne `wc -c <'stdefs.h'`; then
    echo shar: \"'stdefs.h'\" unpacked with wrong size!
fi
# end of 'stdefs.h'
fi
if test -f 'warp.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'warp.c'\"
else
echo shar: Extracting \"'warp.c'\" \(10961 characters\)
sed "s/^X//" >'warp.c' <<'END_OF_FILE'
X
X/*
X * FILE: warp.c
X *   BY: Christopher Lee Fraley ([e-mail])
X * NOTE: Based upon SAIL program by Julius O. Smith (CCRMA, Stanford U)
X * DATE: 17-JUN-88
X */
X
char 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
fail(s)
char *s;
X{
X   printf("warp: %s\n\n", s);     /* Display error message  */
X   exit(1);                       /* Exit, indicating error */
X}
X
fails(s,s2)
char *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
int openData(argc, argv)
int argc;
char *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
closeData()
X{
X   (void) fclose(fin);
X   (void) fclose(fout);
X}
X
X
X
int readData(Data, DataArraySize, Xoff)  /* return: 0 - notDone */
HWORD Data[];                            /*        >0 - index of last sample */
int 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
writeData(Nout, Data)
int Nout;
HWORD Data[];
X{
X   while (--Nout >= 0)                  /* Write Nout samples to ascii file */
X      fprintf(fout, "%d\n", *Data++);
X}
X
X
X
X
check()
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 */
double a,b,c,d,e,f,Total;
int Acc;
X
initWarp()
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
double warpFunction(Time)
UWORD 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
SrcUD(X, Y, Time, Nx, Nwing, LpScl, Imp, ImpD, Interp)
HWORD X[], Y[], Imp[], ImpD[];
UHWORD Nx, Nwing, LpScl;
UWORD *Time;
BOOL 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
int resample(Imp, ImpD, LpScl, Nmult, Nwing, InterpFilt)
HWORD Imp[], ImpD[];
UHWORD LpScl, Nmult, Nwing;
BOOL 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
main(argc,argv)
int argc;
char *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 ([e-mail])\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
END_OF_FILE
if test 10961 -ne `wc -c <'warp.c'`; then
    echo shar: \"'warp.c'\" unpacked with wrong size!
fi
# end of 'warp.c'
fi
echo shar: End of archive 1 \(of 1\).
cp /dev/null ark1isdone
MISSING=""
for I in 1 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have the archive.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0