[xiph-cvs] cvs commit: vorbis-tools/oggenc resample.c resample.h

Michael Smith msmith at xiph.org
Thu May 30 05:14:09 PDT 2002



msmith      02/05/30 05:14:09

  Added:       oggenc   resample.c resample.h
  Log:
  Actually add the resampler.

Revision  Changes    Path
1.1                  vorbis-tools/oggenc/resample.c

Index: resample.c
===================================================================
/* resample.c: see resample.h for interesting stuff */

#include <math.h>
#include <malloc.h>
#include <string.h>
#include <stdarg.h>
#include <assert.h>

#include "resample.h"

<p>static int hcf(int arg1, int arg2)
{
        int mult = 1;

        while (~(arg1 | arg2) & 1)
                arg1 >>= 1, arg2 >>= 1, mult <<= 1;

        while (arg1 > 0)
        {
                if (~(arg1 & arg2) & 1)
                {
                        arg1 >>= (~arg1 & 1);
                        arg2 >>= (~arg2 & 1);
                }
                else if (arg1 < arg2)
                        arg2 = (arg2 - arg1) >> 1;
                else
                        arg1 = (arg1 - arg2) >> 1;
        }

        return arg2 * mult;
}

<p>static void filt_sinc(float *dest, int N, int step, double fc, double gain, int width)
{
        double s = fc / step;
        int mid, x;
        float *endpoint = dest + N,
                *base = dest,
                *origdest = dest;
        
        assert(width <= N);

        if ((N & 1) == 0)
        {
                *dest = 0.0;
                dest += width;
                if (dest >= endpoint)
                        dest = ++base;
                N--;
        }

        mid = N / 2;
        x = -mid;

        while (N--)
        {
                *dest = (x ? sin(x * M_PI * s) / (x * M_PI) * step : fc) * gain;
                x++;
                dest += width;
                if (dest >= endpoint)
                        dest = ++base;
        }
        assert(dest == origdest + width);
}

<p>static double I_zero(double x)
{
        int n = 0;
        double u = 1.0,
                s = 1.0,
                t;

        do
        {
                n += 2;
                t = x / n;
                u *= t * t;
                s += u;
        } while (u > 1e-21 * s);

        return s;
}

<p>static void win_kaiser(float *dest, int N, double alpha, int width)
{
        double I_alpha, midsq;
        int x;
        float *endpoint = dest + N,
                *base = dest,
                *origdest = dest;

        assert(width <= N);

        if ((N & 1) == 0)
        {
                *dest = 0.0;
                dest += width;
                if (dest >= endpoint)
                        dest = ++base;
                N--;
        }

        x = -(N / 2);
        midsq = (double)(x - 1) * (double)(x - 1);
        I_alpha = I_zero(alpha);

        while (N--)
        {
                *dest *= I_zero(alpha * sqrt(1.0 - ((double)x * (double)x) / midsq)) / I_alpha;
                x++;
                dest += width;
                if (dest >= endpoint)
                        dest = ++base;
        }
        assert(dest == origdest + width);
}

<p>int res_init(res_state *state, int channels, int outfreq, int infreq, res_parameter op1, ...)
{
        double beta = 16.0,
                cutoff = 0.80,
                gain = 1.0;
        int taps = 45;

        int factor;

        assert(state);
        assert(channels > 0);
        assert(outfreq > 0);
        assert(infreq > 0);
        assert(taps > 0);

        if (state == NULL || channels <= 0 || outfreq <= 0 || infreq <= 0 || taps <= 0)
                return -1;

        if (op1 != RES_END)
        {
                va_list argp;
                va_start(argp, op1);
                do
                {
                        switch (op1)
                        {
                        case RES_GAIN:
                                gain = va_arg(argp, double);
                                break;

                        case RES_CUTOFF:
                                cutoff = va_arg(argp, double);
                                assert(cutoff > 0.01 && cutoff <= 1.0);
                                break;

                        case RES_TAPS:
                                taps = va_arg(argp, int);
                                assert(taps > 2 && taps < 1000);
                                break;
                                
                        case RES_BETA:
                                beta = va_arg(argp, double);
                                assert(beta > 2.0);
                                break;
                        default:
                                assert("arglist" == "valid");
                                return -1;
                        }
                        op1 = va_arg(argp, res_parameter);
                } while (op1 != RES_END);
                va_end(argp);
        }

        factor = hcf(infreq, outfreq);
        outfreq /= factor;
        infreq /= factor;

        /* adjust to rational values for downsampling */
        if (outfreq < infreq)
        {
                /* push the cutoff frequency down to the output frequency */
                cutoff = cutoff * outfreq / infreq; 

                /* compensate for the sharper roll-off requirement
                 * (this method I found empirically, and don't understand, but it's fast) */
                beta = beta * outfreq * outfreq / (infreq * infreq);
        }

        assert(taps >= (infreq + outfreq - 1) / outfreq);

        if ((state->table = calloc(outfreq * taps, sizeof(float))) == NULL)
                return -1;
        if ((state->pool = calloc(channels * taps, sizeof(SAMPLE))) == NULL)
        {
                free(state->table);
                state->table = NULL;
                return -1;
        }

        state->poolfill = taps / 2 + 1;
        state->channels = channels;
        state->outfreq = outfreq;
        state->infreq = infreq;
        state->taps = taps;
        state->offset = 0;

        filt_sinc(state->table, outfreq * taps, outfreq, cutoff, gain, taps);
        win_kaiser(state->table, outfreq * taps, beta, taps);

        return 0;
}

<p>static SAMPLE sum(float const *scale, int count, SAMPLE const *source, SAMPLE const *trigger, SAMPLE const *reset, int srcstep)
{
        float total = 0.0;

        while (count--)
        {
                total += *source * *scale;

                if (source == trigger)
                        source = reset, srcstep = 1;
                source -= srcstep;
                scale++;
        }

        return total;
}

<p>static int push(res_state const * const state, SAMPLE *pool, int * const poolfill, int * const offset, SAMPLE *dest, int dststep, SAMPLE const *source, int srcstep, size_t srclen)
{
        SAMPLE	* const destbase = dest,
                *poolhead = pool + *poolfill,
                *poolend = pool + state->taps,
                *newpool = pool;
        SAMPLE const *refill, *base, *endpoint;
        int	lencheck;

<p>        assert(state);
        assert(pool);
        assert(poolfill);
        assert(dest);
        assert(source);

        assert(state->poolfill != -1);
        
        lencheck = res_push_check(state, srclen);

        /* fill the pool before diving in */
        while (poolhead < poolend && srclen > 0)
        {
                *poolhead++ = *source;
                source += srcstep;
                srclen--;
        }

        if (srclen <= 0)
                return 0;

        base = source;
        endpoint = source + srclen * srcstep;

        while (source < endpoint)
        {
                *dest = sum(state->table + *offset * state->taps, state->taps, source, base, poolend, srcstep);
                dest += dststep;
                *offset += state->infreq;
                while (*offset >= state->outfreq)
                {
                        *offset -= state->outfreq;
                        source += srcstep;
                }
        }

        assert(dest == destbase + lencheck * dststep);

        /* pretend that source has that underrun data we're not going to get */
        srclen += (source - endpoint) / srcstep;

        /* if we didn't get enough to completely replace the pool, then shift things about a bit */
        if (srclen < state->taps)
        {
                refill = pool + srclen;
                while (refill < poolend)
                        *newpool++ = *refill++;

                refill = source - srclen * srcstep;
        }
        else
                refill = source - state->taps * srcstep;

        /* pull in fresh pool data */
        while (refill < endpoint)
        {
                *newpool++ = *refill;
                refill += srcstep;
        }

        assert(newpool > pool);
        assert(newpool <= poolend);

        *poolfill = newpool - pool;

        return (dest - destbase) / dststep;
}

<p>int res_push_max_input(res_state const * const state, size_t maxoutput)
{
        return maxoutput * state->infreq / state->outfreq;
}

<p>int res_push_check(res_state const * const state, size_t srclen)
{
        if (state->poolfill < state->taps)
                srclen -= state->taps - state->poolfill;

        return (srclen * state->outfreq - state->offset + state->infreq - 1) / state->infreq;
}

<p>int res_push(res_state *state, SAMPLE **dstlist, SAMPLE const **srclist, size_t srclen)
{
        int result = -1, poolfill = -1, offset = -1, i;

        assert(state);
        assert(dstlist);
        assert(srclist);
        assert(state->poolfill >= 0);

        for (i = 0; i < state->channels; i++)
        {
                poolfill = state->poolfill;
                offset = state->offset;
                result = push(state, state->pool + i * state->taps, &poolfill, &offset, dstlist[i], 1, srclist[i], 1, srclen);
        }
        state->poolfill = poolfill;
        state->offset = offset;

        return result;
}

<p>int res_push_interleaved(res_state *state, SAMPLE *dest, SAMPLE const *source, size_t srclen)
{
        int result = -1, poolfill = -1, offset = -1, i;
        
        assert(state);
        assert(dest);
        assert(source);
        assert(state->poolfill >= 0);

        for (i = 0; i < state->channels; i++)
        {
                poolfill = state->poolfill;
                offset = state->offset;
                result = push(state, state->pool + i * state->taps, &poolfill, &offset, dest + i, state->channels, source + i, state->channels, srclen);
        }
        state->poolfill = poolfill;
        state->offset = offset;

        return result;
}

<p>int res_drain(res_state *state, SAMPLE **dstlist)
{
        SAMPLE *tail;
        int result = -1, poolfill = -1, offset = -1, i;

        assert(state);
        assert(dstlist);
        assert(state->poolfill >= 0);

        if ((tail = calloc(state->taps, sizeof(SAMPLE))) == NULL)
                return -1;

        for (i = 0; i < state->channels; i++)
        {
                poolfill = state->poolfill;
                offset = state->offset;
                result = push(state, state->pool + i * state->taps, &poolfill, &offset, dstlist[i], 1, tail, 1, state->taps / 2 - 1);
        }
                
        free(tail);

        state->poolfill = -1;

        return result;
}

<p>int res_drain_interleaved(res_state *state, SAMPLE *dest)
{
        SAMPLE *tail;
        int result = -1, poolfill = -1, offset = -1, i;

        assert(state);
        assert(dest);
        assert(state->poolfill >= 0);

        if ((tail = calloc(state->taps, sizeof(SAMPLE))) == NULL)
                return -1;

        for (i = 0; i < state->channels; i++)
        {
                poolfill = state->poolfill;
                offset = state->offset;
                result = push(state, state->pool + i * state->taps, &poolfill, &offset, dest + i, state->channels, tail, 1, state->taps / 2 - 1);
        }
                
        free(tail);

        state->poolfill = -1;

        return result;
}

<p>void res_clear(res_state *state)
{
        assert(state);
        assert(state->table);
        assert(state->pool);

        free(state->table);
        free(state->pool);
        memset(state, 0, sizeof(*state));
}

<p><p><p>1.1                  vorbis-tools/oggenc/resample.h

Index: resample.h
===================================================================
/* This program is licensed under the GNU Library General Public License,
 * version 2, a copy of which is included with this program (LICENCE.LGPL).
 *   
 * (c) 2002 Simon Hosie <gumboot at clear.net.nz>
 *
 *
 * A resampler
 *
 * reference:
 * 	'Digital Filters', third edition, by R. W. Hamming  ISBN 0-486-65088-X
 *
 * history:
 *	2002-05-31	ready for the world (or some small section thereof)
 *
 *
 * TOOD:
 * 	zero-crossing clipping in coefficient table
 */

#ifndef _RESAMPLE_H_INCLUDED
#define _RESAMPLE_H_INCLUDED

typedef float SAMPLE;

typedef struct
{
        unsigned int channels, infreq, outfreq, taps;
        float *table;
        SAMPLE *pool;

        /* dynamic bits */
        int poolfill;
        int offset;
} res_state;

typedef enum
{
        RES_END,
        RES_GAIN,	/* (double)1.0 */
        RES_CUTOFF,	/* (double)0.80 */ 
        RES_TAPS,	/* (int)45 */
        RES_BETA	/* (double)16.0 */
} res_parameter;

int res_init(res_state *state, int channels, int outfreq, int infreq, res_parameter op1, ...);
/*
 * Configure *state to manage a data stream with the specified parameters.  The
 * string 'params' is currently unspecified, but will configure the parameters
 * of the filter.
 *
 * This function allocates memory, and requires that res_clear() be called when
 * the buffer is no longer needed.
 *
 *
 * All counts/lengths used in the following functions consider only the data in
 * a single channel, and in numbers of samples rather than bytes, even though
 * functionality will be mirrored across as many channels as specified here.
 */

<p>int res_push_max_input(res_state const *state, size_t maxoutput);
/*
 *  Returns the maximum number of input elements that may be provided without
 *  risk of flooding an output buffer of size maxoutput.  maxoutput is
 *  specified in counts of elements, NOT in bytes.
 */

<p>int res_push_check(res_state const *state, size_t srclen);
/*
 * Returns the number of elements that will be returned if the given srclen
 * is used in the next call to res_push().
 */

<p>int res_push(res_state *state, SAMPLE **dstlist, SAMPLE const **srclist, size_t srclen);
int res_push_interleaved(res_state *state, SAMPLE *dest, SAMPLE const *source, size_t srclen);
/*
 * Pushes srclen samples into the front end of the filter, and returns the
 * number of resulting samples.
 *
 * res_push(): srclist and dstlist point to lists of pointers, each of which
 * indicates the beginning of a list of samples.
 *
 * res_push_interleaved(): source and dest point to the beginning of a list of
 * interleaved samples.
 */

<p>int res_drain(res_state *state, SAMPLE **dstlist);
int res_drain_interleaved(res_state *state, SAMPLE *dest);
/*
 * Recover the remaining elements by flushing the internal pool with 0 values,
 * and storing the resulting samples.
 *
 * After either of these functions are called, *state should only re-used in a
 * final call to res_clear().
 */

<p>void res_clear(res_state *state);
/*
 * Free allocated buffers, etc.
 */

#endif

<p><p><p><p>--- >8 ----
List archives:  http://www.xiph.org/archives/
Ogg project homepage: http://www.xiph.org/ogg/
To unsubscribe from this list, send a message to 'cvs-request at xiph.org'
containing only the word 'unsubscribe' in the body.  No subject is needed.
Unsubscribe messages sent to the list will be ignored/filtered.



More information about the commits mailing list