adpcm.c

ADPCM 源码

/*
 * encode.c
 *
 * Usage : encode infile  outfile1 outfile2
 */
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>

/* External function definitions. */
extern void g72x_init_state( struct g72x_state *state_ptr );
extern int g721_encoder( int sl, struct g72x_state *state_ptr, int al, FILE *fpw3 );
extern int g721_decoder( int i, struct g72x_state *state_ptr, int al );
int predictor_zero(struct g72x_state *state_ptr );
int predictor_pole(struct g72x_state *state_ptr );
static int quan( int val, short	*table, int	size );
static int fmult( int an, int srn );
int step_size( struct g72x_state *state_ptr, int al );
int quantize( int d, int y,	short *table, int size );
static int search( int val, short *table, int size );
int reconstruct( int sign, long dqln, int y );
void update( int code_size, int y, long wi, int fi, int dq, int sr, int dqsez, struct g72x_state *state_ptr );
int pack_output( unsigned code, FILE *fpw );
int unpack_input( unsigned char *code, FILE *fpr );

static short seg_end[8] = {0xFF, 0x1FF, 0x3FF, 0x7FF,
				0xFFF, 0x1FFF, 0x3FFF, 0x7FFF};

static short qtab_721[7] = {-124, 80, 178, 246, 300, 349, 400};

/*
 * Maps G.721 code word to reconstructed scale factor normalized log
 * magnitude values.
 */
static long	_dqlntab[16] = {-2048, 4, 135, 213, 273, 323, 373, 425,
				425, 373, 323, 273, 213, 135, 4, -2048};

/* Maps G.721 code word to log of scale factor multiplier. */
static long	_witab[16] = {-12, 18, 41, 64, 112, 198, 355, 1122,
				1122, 355, 198, 112, 64, 41, 18, -12};
/*
 * Maps G.721 code words to a set of values whose long and short
 * term averages are computed and then compared to give an indication
 * how stationary (steady state) the signal is.
 */
static short	_fitab[16] = {0, 0, 0, 0x200, 0x200, 0x200, 0x600, 0xE00,
				0xE00, 0x600, 0x200, 0x200, 0x200, 0, 0, 0};

static short power2[15] = {1, 2, 4, 8, 0x10, 0x20, 0x40, 0x80,
			0x100, 0x200, 0x400, 0x800, 0x1000, 0x2000, 0x4000};

struct g72x_state {
	long yl;	/* Locked or steady state step size multiplier. */
	short yu;	/* Unlocked or non-steady state step size multiplier. */
	short dms;	/* Short term energy estimate. */
	short dml;	/* Long term energy estimate. */
	short ap;	/* Linear weighting coefficient of 'yl' and 'yu'. */

	short a[2];	/* Coefficients of pole portion of prediction filter. */
	short b[6];	/* Coefficients of zero portion of prediction filter. */
	short pk[2];	/*
			 * Signs of previous two samples of a partially
			 * reconstructed signal.
			 */
	short dq[6];	/*
			 * Previous 6 samples of the quantized difference
			 * signal represented in an internal floating point
			 * format.
			 */
	short sr[2];	/*
			 * Previous 2 samples of the quantized difference
			 * signal represented in an internal floating point
			 * format.
			 */
	char td;	/* delayed tone detect, new in 1988 version */
};

/*
 * g72x_init_state()
 *
 * This routine initializes and/or resets the g72x_state structure
 * pointed to by 'state_ptr'.
 * All the initial state values are specified in the CCITT G.721 document.
 */
void
g72x_init_state( state_ptr)
	struct g72x_state *state_ptr;
{
	int		cnta;

	state_ptr->yl = 34816;
	state_ptr->yu = 544;
	state_ptr->dms = 0;
	state_ptr->dml = 0;
	state_ptr->ap = 0;
	for (cnta = 0; cnta < 2; cnta++) {
		state_ptr->a[cnta] = 0;
		state_ptr->pk[cnta] = 0;
		state_ptr->sr[cnta] = 32;
	}
	for (cnta = 0; cnta < 6; cnta++) {
		state_ptr->b[cnta] = 0;
		state_ptr->dq[cnta] = 32;
	}
	state_ptr->td = 0;
}

/*
 * g721_encoder()
 *
 * Encodes the input vale of linear PCM, A-law or u-law data sl and returns
 * the resulting code. -1 is returned for unknown input coding value.
 */
int
g721_encoder( sl, state_ptr, al, fpw3 )
	int		sl;
	struct g72x_state *state_ptr;
	int 		al;
	FILE		*fpw3;
{
	long		wi;
	short		fi;
	short		dq;
	short		sezi, se, sez;		/* ACCUM */
	short		d;			/* SUBTA */
	short		sr;			/* ADDB */
	short		y;			/* MIX */
	short		dqsez;			/* ADDC */
	short		i;
	short		tmp;
//	int			size;

	sl >>= 2;			/* 14-bit dynamic range */
	sezi = predictor_zero(state_ptr);
	sez = sezi >> 1;
	se = (sezi + predictor_pole(state_ptr)) >> 1;	/* estimated signal */

	d = sl - se;				/* estimation difference */
	/* quantize the prediction difference */
	y = step_size(state_ptr, al);		/* quantizer step size */
	i = quantize(d, y, qtab_721, 7);	/* i = ADPCM code */

	dq = reconstruct(i & 8, _dqlntab[i], y);	/* quantized est diff */

	sr = (dq < 0) ? se - (dq & 0x3FFF) : se + dq;	/* reconst. signal */

	dqsez = sr + sez - se;			/* pole prediction diff. */
	update(4, y, _witab[i] << 5, _fitab[i], dq, sr, dqsez, state_ptr);
	return (i);
}

int
g721_decoder( i, state_ptr, al )
	int		i;
	struct g72x_state *state_ptr;
   int 		al;
{
	short		sezi, sei, sez, se;	/* ACCUM */
	short		y;			/* MIX */
	short		sr;			/* ADDB */
	short		dq;
	short		dqsez;

	i &= 0x0f;			/* mask to get proper bits */
	sezi = predictor_zero(state_ptr);
	sez = sezi >> 1;
	sei = sezi + predictor_pole(state_ptr);
	se = sei >> 1;			/* se = estimated signal */

	y = step_size(state_ptr, al);	/* dynamic quantizer step size */

	dq = reconstruct(i & 0x08, _dqlntab[i], y); /* quantized diff. */

	sr = (dq < 0) ? (se - (dq & 0x3FFF)) : se + dq;	/* reconst. signal */

	dqsez = sr - se + sez;			/* pole prediction diff. */

	update(4, y, _witab[i] << 5, _fitab[i], dq, sr, dqsez, state_ptr);

	return (sr << 2);	/* sr was 14-bit dynamic range */
}

/*
 * predictor_zero()
 *
 * computes the estimated signal from 6-zero predictor.
 *
 */
int
predictor_zero( state_ptr )
	struct g72x_state *state_ptr;
{
	int		i;
	int		sezi;

	sezi = fmult(state_ptr->b[0] >> 2, state_ptr->dq[0]);

	for (i = 1; i < 6; i++)			/* ACCUM */
		sezi += fmult(state_ptr->b[i] >> 2, state_ptr->dq[i]);
	return (sezi);
}

/*
 * predictor_pole()
 *
 * computes the estimated signal from 2-pole predictor.
 *
 */
int
predictor_pole( state_ptr )
	struct g72x_state *state_ptr;
{
	return (fmult(state_ptr->a[1] >> 2, state_ptr->sr[1]) +
		fmult(state_ptr->a[0] >> 2, state_ptr->sr[0]));
}

/*
 * quan()
 *
 * quantizes the input val against the table of size short integers.
 * It returns i if table[i - 1] <= val < table[i].
 *
 * Using linear search for simple coding.
 */
static int
quan( val, table, size)
	int		val;
	short		*table;
	int		size;
{
	int		i;

	for (i = 0; i < size; i++)
		if (val < *table++)
			break;
	return (i);
}

/*
 * fmult()
 *
 * returns the integer product of the 14-bit integer "an" and
 * "floating point" representation (4-bit exponent, 6-bit mantessa) "srn".
 */
static int
fmult( an, srn )
	int		an;
	int		srn;
{
	short		anexp;
	short		anmag, anmant;
	short		wanexp, wanmag, wanmant;
	short		retval;

	anmag = (an > 0) ? an : ((-an) & 0x1FFF);

	anexp = quan(anmag, power2, 15) - 6;
	anmant = (anmag == 0) ? 32 :
		(anexp >= 0) ? (long)anmag >> anexp :(long)anmag << -anexp;
	wanexp = anexp + (((long)srn >> 6) & 0xF) - 13;
	wanmant = (anmant * (srn & 0x3f) + 0x30) >> 4;
	retval = (wanexp >= 0) ? (((long)wanmant << wanexp) & 0x7FFF) :
		((long)wanmant >> -wanexp);
	return (((an ^ srn) < 0) ? -retval : retval);
}


/*
 * step_size()
 *
 * computes the quantization step size of the adaptive quantizer.
 *
 */
int
step_size( state_ptr, al )
	struct g72x_state *state_ptr;
   int		al;
{
	int		y;
	int		dif;

	if (state_ptr->ap >= 256)
		return (state_ptr->yu);
	else {
		y = state_ptr->yl >> 6;
		dif = state_ptr->yu - y;
		al = state_ptr->ap >> 2;
		if (dif > 0)
			y += (dif * al) >> 6;
		else if (dif < 0)
			y += (dif * al + 0x3F) >> 6;
		return (y);
	}
}
/*
 * quantize()
 *
 * Given a raw sample, 'd', of the difference signal and a
 * quantization step size scale factor, 'y', this routine returns the
 * ADPCM codeword to which that sample gets quantized.  The step
 * size scale factor division operation is done in the log base 2 domain
 * as a subtraction.
 */
int
quantize( d, y, table, size )
	int		d;	/* Raw difference signal sample */
	int		y;	/* Step size multiplier */
	short		*table;	/* quantization table */
	int		size;	/* table size of short integers */
{
	int		dqm;	/* Magnitude of 'd' */
	short		exp;	/* Integer part of base 2 log of 'd' */
	short		mant;	/* Fractional part of base 2 log */
	short		dl;	/* Log of magnitude of 'd' */
	short		dln;	/* Step size scale factor normalized log */
	int		i;

	/*
	 * LOG
	 *
	 * Compute base 2 log of 'd', and store in 'dl'.
	 */
	dqm = abs(d);
	exp = quan(dqm >> 1, power2, 15);
	mant = (((long)dqm << 7) >> exp) & 0x7F;	/* Fractional portion. */
	dl = ((long)exp << 7) + mant;

	/*
	 * SUBTB
	 *
	 * "Divide" by step size multiplier.
	 */
	dln = dl - (y >> 2);

	/*
	 * QUAN
	 *
	 * Obtain codword i for 'd'.
	 */
	i = quan(dln, table, size);
	if (d < 0)			/* take 1's complement of i */
		return ((size << 1) + 1 - i);
	else if (i == 0)		/* take 1's complement of 0 */
		return ((size << 1) + 1); /* new in 1988 */
	else
		return (i);
}