resample.c

MIDI解码程序(用VC编写)

/*
    TiMidity++ -- MIDI to WAVE converter and player
    Copyright (C) 1999-2002 Masanao Izumo <[e-mail]>
    Copyright (C) 1995 Tuukka Toivonen <[e-mail]>

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    resample.c
*/

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif /* HAVE_CONFIG_H */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#ifndef NO_STRING_H
#include <string.h>
#else
#include <strings.h>
#endif

#include "timidity.h"
#include "common.h"
#include "instrum.h"
#include "playmidi.h"
#include "output.h"
#include "controls.h"
#include "tables.h"
#include "resample.h"
#include "recache.h"


/* for start/end of samples */
static float newt_coeffs[58][58] = {
#include "newton_table.c"
};

int sample_bounds_min, sample_bounds_max; /* min/max bounds for sample data */

/* 4-point interpolation by cubic spline curve. */

static resample_t resample_cspline(sample_t *src, splen_t ofs, resample_rec_t *rec)
{
    int32 ofsi, ofsf, v0, v1, v2, v3, temp;

    ofsi = ofs >> FRACTION_BITS;
    v1 = src[ofsi];
    v2 = src[ofsi + 1];
    if((ofs<rec->loop_start+(1L<<FRACTION_BITS))||
       ((ofs+(2L<<FRACTION_BITS))>rec->loop_end)){
	return (v1 + ((resample_t)((v2 - v1) * (ofs & FRACTION_MASK)) >> FRACTION_BITS));
    } else {
	v0 = src[ofsi - 1];
	v3 = src[ofsi + 2];
	ofsf = ofs & FRACTION_MASK;
	temp = v2;
	v2 = (6 * v2 + ((((5 * v3 - 11 * v2 + 7 * v1 - v0) >> 2) *
			 (ofsf + (1L << FRACTION_BITS)) >> FRACTION_BITS) *
			(ofsf - (1L << FRACTION_BITS)) >> FRACTION_BITS))
	    * ofsf;
	v1 = (((6 * v1+((((5 * v0 - 11 * v1 + 7 * temp - v3) >> 2) *
			 ofsf >> FRACTION_BITS) * (ofsf - (2L << FRACTION_BITS))
			>> FRACTION_BITS)) * ((1L << FRACTION_BITS) - ofsf)) + v2)
	    / (6L << FRACTION_BITS);
	return ((v1 > sample_bounds_max) ? sample_bounds_max :
		((v1 < sample_bounds_min) ? sample_bounds_min : v1));
    }
}


/* 4-point interpolation by Lagrange method.
   Lagrange is now faster than C-spline.  Both have about the same accuracy,
   so choose Lagrange over C-spline, since it is faster.  Technically, it is
   really a 3rd order Newton polynomial (whereas the old Lagrange truely was
   the Lagrange form of the polynomial).  Both Newton and Lagrange forms
   yield the same numerical results, but the Newton form is faster.  Since
   n'th order Newton interpolaiton is resample_newton(), it made sense to
   just keep this labeled as resample_lagrange(), even if it really is the
   Newton form of the polynomial. */

static resample_t resample_lagrange(sample_t *src, splen_t ofs, resample_rec_t *rec)
{
    int32 ofsi, ofsf, v0, v1, v2, v3;

    ofsi = ofs >> FRACTION_BITS;
    v1 = (int32)src[ofsi];
    v2 = (int32)src[ofsi + 1];
    if((ofs<rec->loop_start+(1L<<FRACTION_BITS))||
       ((ofs+(2L<<FRACTION_BITS))>rec->loop_end)) {
	return (v1 + ((resample_t)((v2 - v1) * (ofs & FRACTION_MASK)) >> FRACTION_BITS));
    } else {
	v0 = (int32)src[ofsi - 1];
	v3 = (int32)src[ofsi + 2];
	ofsf = (ofs & FRACTION_MASK) + (1<<FRACTION_BITS);
	v3 += -3*v2 + 3*v1 - v0;
	v3 *= (ofsf - (2<<FRACTION_BITS)) / 6;
	v3 >>= FRACTION_BITS;
	v3 += v2 - v1 - v1 + v0;
	v3 *= (ofsf - (1<<FRACTION_BITS)) >> 1;
	v3 >>= FRACTION_BITS;
	v3 += v1 - v0;
	v3 *= ofsf;
	v3 >>= FRACTION_BITS;
	v3 += v0;
	return ((v3 > sample_bounds_max) ? sample_bounds_max :
		((v3 < sample_bounds_min) ? sample_bounds_min : v3));
    }
}


/* Very fast and accurate table based interpolation.  Better speed and higher
   accuracy than Newton.  This isn't *quite* true Gauss interpolation; it's
   more a slightly modified Gauss interpolation that I accidently stumbled
   upon.  Rather than normalize all x values in the window to be in the range
   [0 to 2*PI], it simply divides them all by 2*PI instead.  I don't know why
   this works, but it does.  Gauss should only work on periodic data with the
   window spanning exactly one period, so it is no surprise that regular Gauss
   interpolation doesn't work too well on general audio data.  But dividing
   the x values by 2*PI magically does.  Any other scaling produces degraded
   results or total garbage.  If anyone can work out the theory behind why
   this works so well (at first glance, it shouldn't ??), please contact me
   (Eric A. Welsh, [e-mail]), as I would really like to have some
   mathematical justification for doing this.  Despite the lack of any sound
   theoretical basis, this method DOES result in highly accurate interpolation
   (or possibly approximaton, not sure yet if it truly interpolates, but it
   looks like it does).  -N 34 is as high as it can go before errors start
   appearing.  But even at -N 34, it is more accurate than Newton at -N 57.
   -N 34 has no problem running in realtime on my system, but -N 25 is the
   default, since that is the optimal compromise between speed and accuracy.
   I strongly recommend using Gauss interpolation.  It is the highest
   quality interpolation option available, and is much faster than using
   Newton polynomials. */

#define DEFAULT_GAUSS_ORDER	25
static float *gauss_table[(1<<FRACTION_BITS)] = {0};	/* don't need doubles */
static int gauss_n = DEFAULT_GAUSS_ORDER;

static resample_t resample_gauss(sample_t *src, splen_t ofs, resample_rec_t *rec)
{
    sample_t *sptr;
    int32 left, right, temp_n;

    left = (ofs>>FRACTION_BITS);
    right = (rec->data_length>>FRACTION_BITS) - left - 1;
    temp_n = (right<<1)-1;
    if (temp_n > (left<<1)+1)
	temp_n = (left<<1)+1;
    if (temp_n < gauss_n) {
	int ii, jj;
	float xd, y;
	if (temp_n <= 0)
	    temp_n = 1;
	xd = ofs & FRACTION_MASK;
	xd /= (1L<<FRACTION_BITS);
	xd += temp_n>>1;
	y = 0;
	sptr = src + (ofs>>FRACTION_BITS) - (temp_n>>1);
	for (ii = temp_n; ii;) {
	    for (jj = 0; jj <= ii; jj++)
		y += sptr[jj] * newt_coeffs[ii][jj];
	    y *= xd - --ii;
	}
	y += *sptr;
	return ((y > sample_bounds_max) ? sample_bounds_max :
		((y < sample_bounds_min) ? sample_bounds_min : y));
    } else {
	float *gptr, *gend;
	float y;
	y = 0;
	sptr = src + left - (gauss_n>>1);
	gptr = gauss_table[ofs&FRACTION_MASK];
	if (gauss_n == DEFAULT_GAUSS_ORDER) {
	    /* expanding the loop for the default case.
	     * this will allow intensive optimization when compiled
	     * with SSE2 capability.
	     */
#define do_gauss  y += *(sptr++) * *(gptr++);
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    do_gauss;
	    y += *sptr * *gptr;
#undef do_gauss
	} else {
	    gend = gptr + gauss_n;
	    do {
		y += *(sptr++) * *(gptr++);
	    } while (gptr <= gend);
	}
	return ((y > sample_bounds_max) ? sample_bounds_max :
		((y < sample_bounds_min) ? sample_bounds_min : y));
    }
}


/* (at least) n+1 point interpolation using Newton polynomials.
   n can be set with a command line option, and
   must be an odd number from 1 to 57 (57 is as high as double precision
   can go without precision errors).  Default n = 11 is good for a 1.533 MHz
   Athlon.  Larger values for n require very fast processors for real time
   playback.  Some points will be interpolated at orders > n to both increase
   accuracy and save CPU. */

static int newt_n = 11;
static int32 newt_old_trunc_x = -1;
static int newt_grow = -1;
static int newt_max = 13;
static double newt_divd[60][60];
static double newt_recip[60] = { 0, 1, 1.0/2, 1.0/3, 1.0/4, 1.0/5, 1.0/6, 1.0/7,
			1.0/8, 1.0/9, 1.0/10, 1.0/11, 1.0/12, 1.0/13, 1.0/14,
			1.0/15, 1.0/16, 1.0/17, 1.0/18, 1.0/19, 1.0/20, 1.0/21,
			1.0/22, 1.0/23, 1.0/24, 1.0/25, 1.0/26, 1.0/27, 1.0/28,
			1.0/29, 1.0/30, 1.0/31, 1.0/32, 1.0/33, 1.0/34, 1.0/35,
			1.0/36, 1.0/37, 1.0/38, 1.0/39, 1.0/40, 1.0/41, 1.0/42,
			1.0/43, 1.0/44, 1.0/45, 1.0/46, 1.0/47, 1.0/48, 1.0/49,
			1.0/50, 1.0/51, 1.0/52, 1.0/53, 1.0/54, 1.0/55, 1.0/56,
			1.0/57, 1.0/58, 1.0/59 };
static sample_t *newt_old_src = NULL;

static resample_t resample_newton(sample_t *src, splen_t ofs, resample_rec_t *rec)
{
    int n_new, n_old;
    int32 v1, v2, diff;
    sample_t *sptr;
    double y, xd;
    int32 left, right, temp_n;
    int ii, jj;

    left = (ofs>>FRACTION_BITS);
    right = (rec->data_length>>FRACTION_BITS)-(ofs>>FRACTION_BITS)-1;
    temp_n = (right<<1)-1;
    if (temp_n <= 0)
	temp_n = 1;
    if (temp_n > (left<<1)+1)
	temp_n = (left<<1)+1;
    if (temp_n < newt_n) {
	xd = ofs & FRACTION_MASK;
	xd /= (1L<<FRACTION_BITS);
	xd += temp_n>>1;
	y = 0;
	sptr = src + (ofs>>FRACTION_BITS) - (temp_n>>1);
	for (ii = temp_n; ii;) {
	    for (jj = 0; jj <= ii; jj++)
		y += sptr[jj] * newt_coeffs[ii][jj];
	    y *= xd - --ii;
	} y += *sptr;
    }else{
	if (newt_grow >= 0 && src == newt_old_src &&
	    (diff = (ofs>>FRACTION_BITS) - newt_old_trunc_x) > 0){
	    n_new = newt_n + ((newt_grow + diff)<<1);
	    if (n_new <= newt_max){
		n_old = newt_n + (newt_grow<<1);
		newt_grow += diff;
		for (v1=(ofs>>FRACTION_BITS)+(n_new>>1)+1,v2=n_new;
		     v2 > n_old; --v1, --v2){
		    newt_divd[0][v2] = src[v1];
		}for (v1 = 1; v1 <= n_new; v1++)
		    for (v2 = n_new; v2 > n_old; --v2)
			newt_divd[v1][v2] = (newt_divd[v1-1][v2] -
					     newt_divd[v1-1][v2-1]) *
			    newt_recip[v1];
	    }else newt_grow = -1;
	}
	if (newt_grow < 0 || src != newt_old_src || diff < 0){
	    newt_grow = 0;
	    for (v1=(ofs>>FRACTION_BITS)-(newt_n>>1),v2=0;
		 v2 <= newt_n; v1++, v2++){
		newt_divd[0][v2] = src[v1];
	    }for (v1 = 1; v1 <= newt_n; v1++)
		for (v2 = newt_n; v2 >= v1; --v2)
		    newt_divd[v1][v2] = (newt_divd[v1-1][v2] -
					 newt_divd[v1-1][v2-1]) *
			newt_recip[v1];
	}
	n_new = newt_n + (newt_grow<<1);
	v2 = n_new;
	y = newt_divd[v2][v2];
	xd = (double)(ofs&FRACTION_MASK) / (1L<<FRACTION_BITS) +
	    (newt_n>>1) + newt_grow;
	for (--v2; v2; --v2){
	    y *= xd - v2;
	    y += newt_divd[v2][v2];
	}y = y*xd + **newt_divd;
	newt_old_src = src;
	newt_old_trunc_x = (ofs>>FRACTION_BITS);
    }
    return ((y > sample_bounds_max) ? sample_bounds_max :
    	    ((y < sample_bounds_min) ? sample_bounds_min : y));
}


/* Simple linear interpolation */

static resample_t resample_linear(sample_t *src, splen_t ofs, resample_rec_t *rec)
{
    int32 v1, v2, ofsi;

    ofsi = ofs >> FRACTION_BITS;
    v1 = src[ofsi];
    v2 = src[ofsi + 1];
#if defined(LOOKUP_HACK) && defined(LOOKUP_INTERPOLATION)
    return (sample_t)(v1 + (iplookup[(((v2 - v1) << 5) & 0x03FE0) |
				     ((ofs & FRACTION_MASK) >> (FRACTION_BITS-5))]));
#else
    return (v1 + ((resample_t)((v2 - v1) * (ofs & FRACTION_MASK)) >> FRACTION_BITS));
#endif
}


/* No interpolation -- Earplugs recommended for maximum listening enjoyment */

static resample_t resample_none(sample_t *src, splen_t ofs, resample_rec_t *rec)
{
    return src[ofs >> FRACTION_BITS];
}


/*
 */

typedef resample_t (*resampler_t)(sample_t*, splen_t, resample_rec_t *);
static resampler_t resamplers[] = {
    resample_cspline,
    resample_lagrange,
    resample_gauss,
    resample_newton,
    resample_linear,
    resample_none
};

#ifdef FIXED_RESAMPLATION
/* don't allow to change the resamplation algorighm.
 * accessing directly to the given function.
 * hope the compiler will optimize the overhead of function calls in this case.
 */
#define cur_resample DEFAULT_RESAMPLATION
#else
static resampler_t cur_resample = DEFAULT_RESAMPLATION;
#endif

#define RESAMPLATION *dest++ = cur_resample(src, ofs, &resrc);

/* exported for recache.c */
resample_t do_resamplation(sample_t *src, splen_t ofs, resample_rec_t *rec)
{
    return cur_resample(src, ofs, rec);
}

/* return the current resampling algorithm */
int get_current_resampler(void)
{
    int i;
    for (i = 0; i < (int)(sizeof(resamplers)/sizeof(resamplers[0])); i++)
	if (resamplers[i] == cur_resample)
	    return i;
    return 0;
}

/* set the current resampling algorithm */
int set_current_resampler(int type)
{
#ifdef FIXED_RESAMPLATION
    return -1;
#else
    if (type < 0 || type > RESAMPLE_NONE)
	return -1;
    cur_resample = resamplers[type];
    return 0;
#endif
}

/* #define FINALINTERP if (ofs < le) *dest++=src[(ofs>>FRACTION_BITS)-1]/2; */
#define FINALINTERP /* Nothing to do after TiMidity++ 2.9.0 */
/* So it isn't interpolation. At least it's final. */

static resample_t resample_buffer[AUDIO_BUFFER_SIZE];
static int resample_buffer_offset;
static resample_t *vib_resample_voice(int, int32 *, int);
static resample_t *normal_resample_voice(int, int32 *, int);

#ifdef PRECALC_LOOPS
#if SAMPLE_LENGTH_BITS == 32 && TIMIDITY_HAVE_INT64
#define PRECALC_LOOP_COUNT(start, end, incr) (int32)(((int64)((end) - (start) + (incr) - 1)) / (incr))
#else
#define PRECALC_LOOP_COUNT(start, end, incr) (int32)(((splen_t)((end) - (start) + (incr) - 1)) / (incr))
#endif
#endif /* PRECALC_LOOPS */

void initialize_gauss_table(int n)
{
    int m, i, k, n_half = (n>>1);
    double ck;
    double x, x_inc, xz;
    double z[35], zsin_[34 + 35], *zsin, xzsin[35];
    float *gptr;

    for (i = 0; i <= n; i++)
    	z[i] = i / (4*M_PI);
    zsin = &zsin_[34];
    for (i = -n; i <= n; i++)
    	zsin[i] = sin(i / (4*M_PI));

    x_inc = 1.0 / (1<<FRACTION_BITS);
    gptr = safe_realloc(gauss_table[0], (n+1)*sizeof(float)*(1<<FRACTION_BITS));
    for (m = 0, x = 0.0; m < (1<<FRACTION_BITS); m++, x += x_inc)
    {
    	xz = (x + n_half) / (4*M_PI);
    	for (i = 0; i <= n; i++)
	    xzsin[i] = sin(xz - z[i]);
    	gauss_table[m] = gptr;

    	for (k = 0; k <= n; k++)
    	{
    	    ck = 1.0;

    	    for (i = 0; i <= n; i++)
    	    {
    	    	if (i == k)
    	    	    continue;
    	
    	    	ck *= xzsin[i] / zsin[k - i];
    	    }
    	    
    	    *gptr++ = ck;
    	}
    }
}

#if 0 /* NOT USED */
/* the was calculated statically in newton_table.c */
static void initialize_newton_coeffs(void)
{
    int i, j, n = 57;
    int sign;

    newt_coeffs[0][0] = 1;
    for (i = 0; i <= n; i++)
    {
    	newt_coeffs[i][0] = 1;
    	newt_coeffs[i][i] = 1;

	if (i > 1)
	{