diff --git a/tools/Makefile b/tools/Makefile index a5f27d1..9ff8d76 100644 --- a/tools/Makefile +++ b/tools/Makefile @@ -24,7 +24,7 @@ aifc_decode_CFLAGS := -O2 # both runs and compiles faster than -O3 aiff_extract_codebook_SOURCES := aiff_extract_codebook.c -tabledesign_SOURCES := sdk-tools/tabledesign/codebook.c sdk-tools/tabledesign/estimate.c sdk-tools/tabledesign/print.c sdk-tools/tabledesign/tabledesign.c +tabledesign_SOURCES := sdk-tools/tabledesign_src/codebook.c sdk-tools/tabledesign_src/estimate.c sdk-tools/tabledesign_src/print.c sdk-tools/tabledesign_src/tabledesign.c tabledesign_CFLAGS := -Wno-uninitialized -laudiofile vadpcm_enc_SOURCES := sdk-tools/adpcm/vadpcm_enc.c sdk-tools/adpcm/vpredictor.c sdk-tools/adpcm/quant.c sdk-tools/adpcm/util.c sdk-tools/adpcm/vencode.c diff --git a/tools/sdk-tools/tabledesign_src/.gitignore b/tools/sdk-tools/tabledesign_src/.gitignore new file mode 100644 index 0000000..04bd425 --- /dev/null +++ b/tools/sdk-tools/tabledesign_src/.gitignore @@ -0,0 +1,9 @@ +*.o +*.s +*.dump +*.aiff +*.aifc +*.table +/tabledesign_irix +/tabledesign_native +/tabledesign_orig diff --git a/tools/sdk-tools/tabledesign_src/Makefile b/tools/sdk-tools/tabledesign_src/Makefile new file mode 100644 index 0000000..a06d70a --- /dev/null +++ b/tools/sdk-tools/tabledesign_src/Makefile @@ -0,0 +1,31 @@ +# Makefile for building tabledesign for either IRIX or natively. +# For an IRIX build, the env variable IRIX_ROOT should point to the root of an +# IRIX filesystem, and QEMU_IRIX should point to the qemu-irix binary. + +IRIX_CC := $(QEMU_IRIX) -silent -L $(IRIX_ROOT) $(IRIX_ROOT)/usr/bin/cc +IRIX_CFLAGS := -fullwarn -Wab,-r4300_mul -Xcpluscomm -mips1 -O2 + +NATIVE_CC := gcc +NATIVE_CFLAGS := -Wall -Wno-uninitialized -O2 + +LDFLAGS := -lm -laudiofile + +default: native +all: irix native + +irix: tabledesign_irix +native: tabledesign_native + +clean: + $(RM) *.o tabledesign_irix tabledesign_native + +%.o: %.c + $(IRIX_CC) -c $(IRIX_CFLAGS) $< -o $@ + +tabledesign_irix: tabledesign.o codebook.o estimate.o print.o + $(IRIX_CC) $^ -o $@ $(LDFLAGS) + +tabledesign_native: tabledesign.c codebook.c estimate.c print.c + $(NATIVE_CC) $(NATIVE_CFLAGS) $^ -o $@ $(LDFLAGS) + +.PHONY: default all irix native clean diff --git a/tools/sdk-tools/tabledesign_src/codebook.c b/tools/sdk-tools/tabledesign_src/codebook.c new file mode 100644 index 0000000..579a168 --- /dev/null +++ b/tools/sdk-tools/tabledesign_src/codebook.c @@ -0,0 +1,104 @@ +#include +#include "tabledesign.h" + +void split(double **table, double *delta, int order, int npredictors, double scale) +{ + int i, j; + + for (i = 0; i < npredictors; i++) + { + for (j = 0; j <= order; j++) + { + table[i + npredictors][j] = table[i][j] + delta[j] * scale; + } + } +} + +void refine(double **table, int order, int npredictors, double **data, int dataSize, int refineIters, UNUSED double unused) +{ + int iter; // spD8 + double **rsums; + int *counts; // spD0 + double *temp_s7; + double dist; + double dummy; // spC0 + double bestValue; + int bestIndex; + int i, j; + + rsums = malloc(npredictors * sizeof(double*)); + for (i = 0; i < npredictors; i++) + { + rsums[i] = malloc((order + 1) * sizeof(double)); + } + + counts = malloc(npredictors * sizeof(int)); + temp_s7 = malloc((order + 1) * sizeof(double)); + + for (iter = 0; iter < refineIters; iter++) + { + for (i = 0; i < npredictors; i++) + { + counts[i] = 0; + for (j = 0; j <= order; j++) + { + rsums[i][j] = 0.0; + } + } + + for (i = 0; i < dataSize; i++) + { + bestValue = 1e30; + bestIndex = 0; + + for (j = 0; j < npredictors; j++) + { + dist = model_dist(table[j], data[i], order); + if (dist < bestValue) + { + bestValue = dist; + bestIndex = j; + } + } + + counts[bestIndex]++; + rfroma(data[i], order, temp_s7); + for (j = 0; j <= order; j++) + { + rsums[bestIndex][j] += temp_s7[j]; + } + } + + for (i = 0; i < npredictors; i++) + { + if (counts[i] > 0) + { + for (j = 0; j <= order; j++) + { + rsums[i][j] /= counts[i]; + } + } + } + + for (i = 0; i < npredictors; i++) + { + durbin(rsums[i], order, temp_s7, table[i], &dummy); + + for (j = 1; j <= order; j++) + { + if (temp_s7[j] >= 1.0) temp_s7[j] = 0.9999999999; + if (temp_s7[j] <= -1.0) temp_s7[j] = -0.9999999999; + } + + afromk(temp_s7, table[i], order); + } + } + + free(counts); + for (i = 0; i < npredictors; i++) + { + free(rsums[i]); + } + free(rsums); + free(temp_s7); +} diff --git a/tools/sdk-tools/tabledesign_src/estimate.c b/tools/sdk-tools/tabledesign_src/estimate.c new file mode 100644 index 0000000..da2b4e8 --- /dev/null +++ b/tools/sdk-tools/tabledesign_src/estimate.c @@ -0,0 +1,342 @@ +#include +#include +#include "tabledesign.h" + +/** + * Computes the autocorrelation of a vector. More precisely, it computes the + * dot products of vec[i:] and vec[:-i] for i in [0, k). Unused. + * + * See https://en.wikipedia.org/wiki/Autocorrelation. + */ +void acf(double *vec, int n, double *out, int k) +{ + int i, j; + double sum; + for (i = 0; i < k; i++) + { + sum = 0.0; + for (j = 0; j < n - i; j++) + { + sum += vec[j + i] * vec[j]; + } + out[i] = sum; + } +} + +// https://en.wikipedia.org/wiki/Durbin%E2%80%93Watson_statistic ? +// "detects the presence of autocorrelation at lag 1 in the residuals (prediction errors)" +int durbin(double *arg0, int n, double *arg2, double *arg3, double *outSomething) +{ + int i, j; + double sum, div; + int ret; + + arg3[0] = 1.0; + div = arg0[0]; + ret = 0; + + for (i = 1; i <= n; i++) + { + sum = 0.0; + for (j = 1; j <= i-1; j++) + { + sum += arg3[j] * arg0[i - j]; + } + + arg3[i] = (div > 0.0 ? -(arg0[i] + sum) / div : 0.0); + arg2[i] = arg3[i]; + + if (fabs(arg2[i]) > 1.0) + { + ret++; + } + + for (j = 1; j < i; j++) + { + arg3[j] += arg3[i - j] * arg3[i]; + } + + div *= 1.0 - arg3[i] * arg3[i]; + } + *outSomething = div; + return ret; +} + +void afromk(double *in, double *out, int n) +{ + int i, j; + out[0] = 1.0; + for (i = 1; i <= n; i++) + { + out[i] = in[i]; + for (j = 1; j <= i - 1; j++) + { + out[j] += out[i - j] * out[i]; + } + } +} + +int kfroma(double *in, double *out, int n) +{ + int i, j; + double div; + double temp; + double *next; + int ret; + + ret = 0; + next = malloc((n + 1) * sizeof(double)); + + out[n] = in[n]; + for (i = n - 1; i >= 1; i--) + { + for (j = 0; j <= i; j++) + { + temp = out[i + 1]; + div = 1.0 - (temp * temp); + if (div == 0.0) + { + free(next); + return 1; + } + next[j] = (in[j] - in[i + 1 - j] * temp) / div; + } + + for (j = 0; j <= i; j++) + { + in[j] = next[j]; + } + + out[i] = next[i]; + if (fabs(out[i]) > 1.0) + { + ret++; + } + } + + free(next); + return ret; +} + +void rfroma(double *arg0, int n, double *arg2) +{ + int i, j; + double **mat; + double div; + + mat = malloc((n + 1) * sizeof(double*)); + mat[n] = malloc((n + 1) * sizeof(double)); + mat[n][0] = 1.0; + for (i = 1; i <= n; i++) + { + mat[n][i] = -arg0[i]; + } + + for (i = n; i >= 1; i--) + { + mat[i - 1] = malloc(i * sizeof(double)); + div = 1.0 - mat[i][i] * mat[i][i]; + for (j = 1; j <= i - 1; j++) + { + mat[i - 1][j] = (mat[i][i - j] * mat[i][i] + mat[i][j]) / div; + } + } + + arg2[0] = 1.0; + for (i = 1; i <= n; i++) + { + arg2[i] = 0.0; + for (j = 1; j <= i; j++) + { + arg2[i] += mat[i][j] * arg2[i - j]; + } + } + + free(mat[n]); + for (i = n; i > 0; i--) + { + free(mat[i - 1]); + } + free(mat); +} + +double model_dist(double *arg0, double *arg1, int n) +{ + double *sp3C; + double *sp38; + double ret; + int i, j; + + sp3C = malloc((n + 1) * sizeof(double)); + sp38 = malloc((n + 1) * sizeof(double)); + rfroma(arg1, n, sp3C); + + for (i = 0; i <= n; i++) + { + sp38[i] = 0.0; + for (j = 0; j <= n - i; j++) + { + sp38[i] += arg0[j] * arg0[i + j]; + } + } + + ret = sp38[0] * sp3C[0]; + for (i = 1; i <= n; i++) + { + ret += 2 * sp3C[i] * sp38[i]; + } + + free(sp3C); + free(sp38); + return ret; +} + +// compute autocorrelation matrix? +void acmat(short *in, int n, int m, double **out) +{ + int i, j, k; + for (i = 1; i <= n; i++) + { + for (j = 1; j <= n; j++) + { + out[i][j] = 0.0; + for (k = 0; k < m; k++) + { + out[i][j] += in[k - i] * in[k - j]; + } + } + } +} + +// compute autocorrelation vector? +void acvect(short *in, int n, int m, double *out) +{ + int i, j; + for (i = 0; i <= n; i++) + { + out[i] = 0.0; + for (j = 0; j < m; j++) + { + out[i] -= in[j - i] * in[j]; + } + } +} + +/** + * Replaces a real n-by-n matrix "a" with the LU decomposition of a row-wise + * permutation of itself. + * + * Input parameters: + * a: The matrix which is operated on. 1-indexed; it should be of size + * (n+1) x (n+1), and row/column index 0 is not used. + * n: The size of the matrix. + * + * Output parameters: + * indx: The row permutation performed. 1-indexed; it should be of size n+1, + * and index 0 is not used. + * d: the determinant of the permutation matrix. + * + * Returns 1 to indicate failure if the matrix is singular or has zeroes on the + * diagonal, 0 on success. + * + * Derived from ludcmp in "Numerical Recipes in C: The Art of Scientific Computing", + * with modified error handling. + */ +int lud(double **a, int n, int *indx, int *d) +{ + int i,imax,j,k; + double big,dum,sum,temp; + double min,max; + double *vv; + + vv = malloc((n + 1) * sizeof(double)); + *d=1; + for (i=1;i<=n;i++) { + big=0.0; + for (j=1;j<=n;j++) + if ((temp=fabs(a[i][j])) > big) big=temp; + if (big == 0.0) return 1; + vv[i]=1.0/big; + } + for (j=1;j<=n;j++) { + for (i=1;i= big) { + big=dum; + imax=i; + } + } + if (j != imax) { + for (k=1;k<=n;k++) { + dum=a[imax][k]; + a[imax][k]=a[j][k]; + a[j][k]=dum; + } + *d = -(*d); + vv[imax]=vv[j]; + } + indx[j]=imax; + if (a[j][j] == 0.0) return 1; + if (j != n) { + dum=1.0/(a[j][j]); + for (i=j+1;i<=n;i++) a[i][j] *= dum; + } + } + free(vv); + + min = 1e10; + max = 0.0; + for (i = 1; i <= n; i++) + { + temp = fabs(a[i][i]); + if (temp < min) min = temp; + if (temp > max) max = temp; + } + return min / max < 1e-10 ? 1 : 0; +} + +/** + * Solves the set of n linear equations Ax = b, using LU decomposition + * back-substitution. + * + * Input parameters: + * a: The LU decomposition of a matrix, created by "lud". + * n: The size of the matrix. + * indx: Row permutation vector, created by "lud". + * b: The vector b in the equation. 1-indexed; is should be of size n+1, and + * index 0 is not used. + * + * Output parameters: + * b: The output vector x. 1-indexed. + * + * From "Numerical Recipes in C: The Art of Scientific Computing". + */ +void lubksb(double **a, int n, int *indx, double *b) +{ + int i,ii=0,ip,j; + double sum; + + for (i=1;i<=n;i++) { + ip=indx[i]; + sum=b[ip]; + b[ip]=b[i]; + if (ii) + for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; + else if (sum) ii=i; + b[i]=sum; + } + for (i=n;i>=1;i--) { + sum=b[i]; + for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j]; + b[i]=sum/a[i][i]; + } +} diff --git a/tools/sdk-tools/tabledesign_src/print.c b/tools/sdk-tools/tabledesign_src/print.c new file mode 100644 index 0000000..80fe9e0 --- /dev/null +++ b/tools/sdk-tools/tabledesign_src/print.c @@ -0,0 +1,89 @@ +#include +#include +#include "tabledesign.h" + +int print_entry(FILE *out, double *row, int order) +{ + double **table; + double fval; + int ival; + int i, j, k; + int overflows; + + table = malloc(8 * sizeof(double*)); + + for (i = 0; i < 8; i++) + { + table[i] = malloc(order * sizeof(double)); + } + + for (i = 0; i < order; i++) + { + for (j = 0; j < i; j++) + { + table[i][j] = 0.0; + } + + for (j = i; j < order; j++) + { + table[i][j] = -row[order - j + i]; + } + } + + for (i = order; i < 8; i++) + { + for (j = 0; j < order; j++) + { + table[i][j] = 0.0; + } + } + + for (i = 1; i < 8; i++) + { + for (j = 1; j <= order; j++) + { + if (i - j >= 0) + { + for (k = 0; k < order; k++) + { + table[i][k] -= row[j] * table[i - j][k]; + } + } + } + } + + overflows = 0; + for (i = 0; i < order; i++) + { + for (j = 0; j < 8; j++) + { + fval = table[j][i] * 2048.0; + if (fval < 0.0) + { + ival = (int) (fval - 0.5); + if (ival < -0x8000) + { + overflows++; + } + } + else + { + ival = (int) (fval + 0.5); + if (ival >= 0x8000) + { + overflows++; + } + } + fprintf(out, "%5d ", ival); + } + + fprintf(out, "\n"); + } + + for (i = 0; i < 8; i++) + { + free(table[i]); + } + free(table); + return overflows; +} diff --git a/tools/sdk-tools/tabledesign_src/tabledesign.c b/tools/sdk-tools/tabledesign_src/tabledesign.c new file mode 100644 index 0000000..f85757c --- /dev/null +++ b/tools/sdk-tools/tabledesign_src/tabledesign.c @@ -0,0 +1,262 @@ +#include +#include +#include +#include +#include +#include "tabledesign.h" + +#ifdef __sgi + +typedef long SampleFormat; + +#define MODE_READ "r" + +#else + +// The modern implementation of SGI's audiofile library which is in Ubuntu +// (https://github.com/mpruett/audiofile/) has renamed some of the functions, +// and changed some data types. + +typedef int SampleFormat; +#define AFopenfile afOpenFile +#define AFgetchannels afGetChannels +#define AFgettrackids afGetTrackIDs +#define AFgetsampfmt afGetSampleFormat +#define AFgetframecnt afGetFrameCount +#define AFgetrate afGetRate +#define AFreadframes afReadFrames + +#define MODE_READ "rb" + +#endif + +char usage[80] = "[-o order -s bits -t thresh -i refine_iter -f frame_size] aifcfile"; + +int main(int argc, char **argv) +{ + const char *programName; // sp118 + double thresh; // sp110 + int order; // sp10C + int bits; // sp108 + int refineIters; // sp104 + int frameSize; // sp100 + UNUSED int rate; + int frameCount; + int opt; + double *spF4; + double dummy; // spE8 + double **mat; // spE4 + double **data; // spD0 + double *splitDelta; // spCC + int j; // spC0 + int permDet; + int curBits; // spB8 + int npredictors; // spB4 + int *perm; // spB0 + int numOverflows; // spAC + SampleFormat sampleFormat; // sp90 + SampleFormat sampleWidth; // sp8C + AFfilehandle afFile; // sp88 + int channels; + int tracks; + double *vec; // s2 + double **temp_s1; + short *temp_s3; + int i; + int dataSize; // s4 + + order = 2; + bits = 2; + refineIters = 2; + frameSize = 16; + numOverflows = 0; + programName = argv[0]; + thresh = 10.0; + + if (argc < 2) + { + fprintf(stderr, "%s %s\n", argv[0], usage); + exit(1); + } + + while ((opt = getopt(argc, argv, "o:s:t:i:f:")) != -1) + { + switch (opt) + { + case 'o': + if (sscanf(optarg, "%d", &order) != 1) + order = 2; + break; + case 's': + if (sscanf(optarg, "%d", &bits) != 1) + bits = 2; + break; + case 'f': + if (sscanf(optarg, "%d", &frameSize) != 1) + frameSize = 16; + break; + case 'i': + if (sscanf(optarg, "%d", &refineIters) != 1) + refineIters = 2; + break; + case 't': + if (sscanf(optarg, "%lf", &thresh) != 1) + thresh = 10.0; + break; + } + } + + argv = &argv[optind - 1]; + + afFile = AFopenfile(argv[1], MODE_READ, NULL); + if (afFile == NULL) + { + fprintf(stderr, + "%s: input AIFC file [%s] could not be opened.\n", + programName, argv[1]); + exit(1); + } + + channels = AFgetchannels(afFile, AF_DEFAULT_TRACK); + if (channels != 1) + { + fprintf(stderr, + "%s: file [%s] contains %d channels, only 1 channel supported.\n", + programName, argv[1], channels); + exit(1); + } + + tracks = AFgettrackids(afFile, NULL); + if (tracks != 1) + { + fprintf(stderr, + "%s: file [%s] contains %d tracks, only 1 track supported.\n", + programName, argv[1], tracks); + exit(1); + } + + AFgetsampfmt(afFile, AF_DEFAULT_TRACK, &sampleFormat, &sampleWidth); + if (sampleWidth != 16) + { + fprintf(stderr, + "%s: file [%s] contains %d bit samples, only 16 bit samples supported.\n", + programName, argv[1], (int)sampleWidth); + exit(1); + } + + temp_s1 = malloc((1 << bits) * sizeof(double*)); + for (i = 0; i < (1 << bits); i++) + { + temp_s1[i] = malloc((order + 1) * sizeof(double)); + } + + splitDelta = malloc((order + 1) * sizeof(double)); + temp_s3 = malloc(frameSize * 2 * sizeof(short)); + for (i = 0; i < frameSize * 2; i++) + { + temp_s3[i] = 0; + } + + vec = malloc((order + 1) * sizeof(double)); + spF4 = malloc((order + 1) * sizeof(double)); + mat = malloc((order + 1) * sizeof(double*)); + for (i = 0; i <= order; i++) + { + mat[i] = malloc((order + 1) * sizeof(double)); + } + + perm = malloc((order + 1) * sizeof(int)); + frameCount = AFgetframecnt(afFile, AF_DEFAULT_TRACK); + rate = AFgetrate(afFile, AF_DEFAULT_TRACK); + data = malloc(frameCount * sizeof(double*)); + dataSize = 0; + + while (AFreadframes(afFile, AF_DEFAULT_TRACK, temp_s3 + frameSize, frameSize) == frameSize) + { + acvect(temp_s3 + frameSize, order, frameSize, vec); + if (fabs(vec[0]) > thresh) + { + acmat(temp_s3 + frameSize, order, frameSize, mat); + if (lud(mat, order, perm, &permDet) == 0) + { + lubksb(mat, order, perm, vec); + vec[0] = 1.0; + if (kfroma(vec, spF4, order) == 0) + { + data[dataSize] = malloc((order + 1) * sizeof(double)); + data[dataSize][0] = 1.0; + + for (i = 1; i <= order; i++) + { + if (spF4[i] >= 1.0) spF4[i] = 0.9999999999; + if (spF4[i] <= -1.0) spF4[i] = -0.9999999999; + } + + afromk(spF4, data[dataSize], order); + dataSize++; + } + } + } + + for (i = 0; i < frameSize; i++) + { + temp_s3[i] = temp_s3[i + frameSize]; + } + } + + vec[0] = 1.0; + for (j = 1; j <= order; j++) + { + vec[j] = 0.0; + } + + for (i = 0; i < dataSize; i++) + { + rfroma(data[i], order, temp_s1[0]); + for (j = 1; j <= order; j++) + { + vec[j] += temp_s1[0][j]; + } + } + + for (j = 1; j <= order; j++) + { + vec[j] /= dataSize; + } + + durbin(vec, order, spF4, temp_s1[0], &dummy); + + for (j = 1; j <= order; j++) + { + if (spF4[j] >= 1.0) spF4[j] = 0.9999999999; + if (spF4[j] <= -1.0) spF4[j] = -0.9999999999; + } + + afromk(spF4, temp_s1[0], order); + curBits = 0; + while (curBits < bits) + { + for (i = 0; i <= order; i++) + { + splitDelta[i] = 0.0; + } + splitDelta[order - 1] = -1.0; + split(temp_s1, splitDelta, order, 1 << curBits, 0.01); + curBits++; + refine(temp_s1, order, 1 << curBits, data, dataSize, refineIters, 0.0); + } + + npredictors = 1 << curBits; + fprintf(stdout, "%d\n%d\n", order, npredictors); + + for (i = 0; i < npredictors; i++) + { + numOverflows += print_entry(stdout, temp_s1[i], order); + } + + if (numOverflows > 0) + { + fprintf(stderr, "There was overflow - check the table\n"); + } + return 0; +} diff --git a/tools/sdk-tools/tabledesign_src/tabledesign.h b/tools/sdk-tools/tabledesign_src/tabledesign.h new file mode 100644 index 0000000..89c7a15 --- /dev/null +++ b/tools/sdk-tools/tabledesign_src/tabledesign.h @@ -0,0 +1,30 @@ +#ifndef TABLEDESIGN_H +#define TABLEDESIGN_H + +#include + +#ifdef __GNUC__ +#define UNUSED __attribute__((unused)) +#else +#define UNUSED +#endif + +// estimate.c +int durbin(double *thing, int n, double *thing2, double *thing3, double *outSomething); +void afromk(double *in, double *out, int n); +int kfroma(double *in, double *out, int n); +void rfroma(double *dataRow, int n, double *thing3); +double model_dist(double *first, double *second, int n); +void acmat(short *in, int n, int m, double **mat); +void acvect(short *in, int n, int m, double *vec); +int lud(double **a, int n, int *indx, int *d); +void lubksb(double **a, int n, int *indx, double *b); + +// codebook.c +void split(double **table, double *delta, int order, int npredictors, double scale); +void refine(double **table, int order, int npredictors, double **data, int dataSize, int refineIters, double unused); + +// print.c +int print_entry(FILE *out, double *row, int order); + +#endif