/*
 * Copyright (c) 1991 Regents of the University of California.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms are permitted
 * provided that the above copyright notice and this paragraph are
 * duplicated in all such forms and that any documentation,
 * advertising materials, and other materials related to such
 * distribution and use acknowledge that the software was developed
 * by the University of California, Lawrence Berkeley Laboratory,
 * Berkeley, CA.  The name of the University may not be used to
 * endorse or promote products derived from this software without
 * specific prior written permission.
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
 * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 */
static const char rcsid[] =
    "@(#) $Header: mktab.c,v 1.5 96/03/18 00:52:01 van Exp $ (LBL)";

#include "config.h"
#include <stdio.h>
#include <math.h>
#include "mulaw.h"

#define BIAS 0x84
#define CLIP 32635

u_int8_t
st_linear_to_ulaw(register int linear)
{
	static int exp_lut[256] = {0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,
				   4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
				   5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
				   5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
				   6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
				   6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
				   6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
				   6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
				   7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
				   7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
				   7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
				   7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
				   7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
				   7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
				   7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
				   7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7};
	register int sign, exponent, mantissa;
	register u_int8_t ulawbyte;
	
	/* Get the value into sign-magnitude. */
	sign = (linear >> 8) & 0x80;
	if (sign != 0) {
		/* sign extend the value, just in case */
		linear |= 0xffff0000;
		linear = -linear;
	}
	if (linear > CLIP)
		linear = CLIP;
	
	/* Convert from 16 bit linear to ulaw. */
	linear += BIAS;
	exponent = exp_lut[(linear >> 7) & 0xFF];
	mantissa = (linear >> (exponent + 3)) & 0x0F;
	ulawbyte = ~(sign | (exponent << 4) | mantissa);
	return ulawbyte;
}

int
st_ulaw_to_linear(u_int8_t ulawbyte)
{
	static int exp_lut[8] = { 0, 132, 396, 924, 1980, 4092, 8316, 16764 };
	int sign, exponent, mantissa, sample;
	
	ulawbyte = ~ulawbyte;
	sign = (ulawbyte & 0x80);
	exponent = (ulawbyte >> 4) & 0x07;
	mantissa = ulawbyte & 0x0F;
	sample = exp_lut[exponent] + (mantissa << (exponent + 3)) +
		 (1 << (exponent + 2));
	if (sign != 0)
		sample = -sample;
	/*
	 * ulaw maps zero to full-scale negative.  This is a very poor
	 * choice given that most i/o failure modes stick zero in the
	 * buffer.  So, to save our ears, we map ulaw zero to linear zero.
	 */
	return (ulawbyte == 0xff? 0 : sample);
}

void
tab_tolin(register FILE *f)
{
	register int i, j;

	fprintf(f, "const int16_t mulawtolin[256] = {\n");
	for (i = 0; i < 256; ) {
		putc('\t', f);
		for (j = i + 8; i < j; ++i) 
			fprintf(f, "%d, ", st_ulaw_to_linear(i));
		putc('\n', f);
	}
	fprintf(f, "};\n");
}

void
tab_tomu(register FILE *f)
{
	register int i, j;

	fprintf(f, "const u_int8_t lintomulaw[65536] = {\n");
	for (i = 0; i < 65536; ) {
		putc('\t', f);
		for (j = i + 8; i < j; ++i) 
			fprintf(f, "%d, ", st_linear_to_ulaw(i));
		putc('\n', f);
	}
	fprintf(f, "};\n");
}

/*
 * build a linear to mulaw lookup table extended by 1
 * bit to handle clipping in the table lookup rather
 * than testing each sample.
 */
void
tab_tomuX(register FILE *f)
{
	register int i, j;

	fprintf(f, "const u_int8_t lintomulawX[65536*2] = {\n");
	/* positive half */
	for (i = 0; i < 32768; ) {
		putc('\t', f);
		for (j = i + 8; i < j; ++i) 
			fprintf(f, "%d, ", st_linear_to_ulaw(i));
		putc('\n', f);
	}
	/* positive overflow (clipped) */
	for (i = 32768; i < 65536; ) {
		putc('\t', f);
		for (j = i + 8; i < j; ++i) 
			fprintf(f, "%d, ", 0x80);
		putc('\n', f);
	}
	/* negative overflow (clipped) */
	for (i = 65536; i < 98304; ) {
		putc('\t', f);
		for (j = i + 8; i < j; ++i) 
			fprintf(f, "%d, ", 0x01);
		putc('\n', f);
	}
	/* negative half */
	for (i = 98304; i < 131072; ) {
		putc('\t', f);
		for (j = i + 8; i < j; ++i) 
			fprintf(f, "%d, ", st_linear_to_ulaw(i));
		putc('\n', f);
	}
	fprintf(f, "};\n");
}

void
prtone(register FILE *f, register int u)
{
	if (u < -CLIP)
		u = -CLIP;
	else if (u > CLIP)
		u = CLIP;
	fprintf(f, "%d, ", st_linear_to_ulaw(u));
}

void
tab_sum(register FILE *f)
{
	register int i, j, k, u;

	fprintf(f, "const u_int8_t mulawsum[256*256] = {\n");
	for (i = 0; i < 256; ++i) {
		fprintf(f, "/* %d */\n", i);
		u = st_ulaw_to_linear(i);
		for (j = 0; j < 256; ) {
			putc('\t', f);
			for (k = j + 8; j < k; ++j)
				prtone(f, u + st_ulaw_to_linear(j));
			putc('\n', f);
		}
	}
	fprintf(f, "};\n");
}

void
prt8(register FILE *f, register int u, register int* us)
{
	register int k;

	putc('\t', f);
	for (k = 0; k < 8; ++k)
		prtone(f, u + us[k]);
	putc('\n', f);
}

void
prtup64(register FILE *f, register int u, register int uj)
{
	register int k;
	int us[64], uabs, uk;

	/*
	 * compute the 64 values that would precede uj in a mulaw
	 * table (assuming a gain of 17/16 which is roughly right
	 * for large uj).  Do them in reverse order since mulaw
	 * values are complemented.
	 */
	uabs = (uj>=0? uj : -uj) << 4;
	for (k = 64; --k >= 0; ) {
		uk = uabs >> 4;
		us[k] = uj>=0? uk : -uk;
		uabs += uk;
	}
	for (k = 0; k < 64; k += 8)
		prt8(f, u, us + k);
}

void
prtlow64(register FILE *f, register int u, register int uj)
{
	register int k;
	int us[64], uabs, uk;

	/*
	 * compute the 64 values that would follow uj in a mulaw
	 * table.
	 */
	uabs = 1 << 4;
	for (k = 48; --k >= 0; ) {
		uk = uabs >> 4;
		us[k] = uj>=0? uk : -uk;
		uabs += uk;
	}
	for (k = 64; --k >= 48; ) {
		us[k] = 0;
	}
	for (k = 0; k < 64; k += 8)
		prt8(f, u, us + k);
}

void
tab_mugain(register FILE *f)
{
	register int i, j, k;
#define M 16.
	double g = 1.;
	double gi = pow(M, (double)1. / 32.);

	fprintf(f, "const u_int8_t mugaintab[64*256] = {\n");
	for (i = 0; i < 64; ++i) {
		fprintf(f, "/* step %d (gain %g) */\n", i - 32, g / M);

		for (j = 0; j < 256; ) {
			putc('\t', f);
			for (k = j + 8; j < k; ++j) {
				double v = g / M * (double)st_ulaw_to_linear(j);
				prtone(f, (int)v);
			}
			putc('\n', f);
		}
		g *= gi;
	}
	fprintf(f, "};\n");
#undef M
}

/*
 * Output a table that maps mulaw samples to power levels, where
 * the power level is given by p = log_2(x^2) (i.e., not really db).
 * p is expressed in 32-bit fixed point, with 19 bits of fraction.
 * Since ulaw samples range over 16 bits (linear), we need 32 bits
 * to represent log_2(x^2).  To compute the power in a frame, we add
 * 160 of these samples, so we need to be able to represent 32 * 160 = 0x1400,
 * with the encoding.  Thus, we need 13 bits for the integer part,
 * which gives 19 bits of fraction.  This encoding will allow the
 * frame size to go to 256.  Note that all the values are positive
 * since the smallest linear value is 4 (except for 0, which we will
 * just encode as 0).
 */
#ifdef notdef
tab_mudb(f)
	register FILE *f;
{
	register int i, j;
	double log2 = log((double)2);

	fprintf(f, "const unsigned long mulawtodb[256] = {\n");
	for (i = 0; i < 256; ) {
		putc('\t', f);
		for (j = i + 4; i < j; ++i) {
			u_long s;
			double d = (double)st_ulaw_to_linear(i);
			if (d != 0) {
				d = d * d;
				d = log(d) / log2;
				d *= (double)(1 << MUDB_FRACBITS);
			}
			fprintf(f, "0x%lx, ", (u_long)irint(d));
		}
		putc('\n', f);
	}
	fprintf(f, "};\n");
}

tab_x(f)
	register FILE *f;
{
	register int i;

	for (i = 0; i <= CLIP; ++i) {
		int u = st_linear_to_ulaw(i);
		printf("%d\n", ~u & 0x7f);
	}
}

tab_xx(f)
	register FILE *f;
{
	register int i;

	for (i = 0; i <= CLIP; ++i) {
		int u = st_linear_to_ulaw(i);
		printf("%d\n", st_ulaw_to_linear(u));
	}
}
#endif

int
main(int argc, char* argv[])
{
	if (argc != 2)
		exit(1);
	fprintf(stdout, "#include \"config.h\"\n");
	if (strcmp(argv[1], "-mulaw") == 0) {
		tab_tolin(stdout);
		tab_tomu(stdout);
	} else if (strcmp(argv[1], "-sum") == 0)
		tab_sum(stdout);
	else if (strcmp(argv[1], "-scale") == 0)
		tab_mugain(stdout);
	else if (strcmp(argv[1], "-muX") == 0)
		tab_tomuX(stdout);
#ifdef notdef
	else if (strcmp(argv[1], "-db") == 0)
		tab_mudb(stdout);
	else if (strcmp(argv[1], "-x") == 0)
		tab_x(stdout);
	else if (strcmp(argv[1], "-xx") == 0)
		tab_xx(stdout);
#endif
	else
		exit(1);
	return 0;
}