fix tabledesign

master
n64 2019-08-25 00:58:46 -04:00
parent 89e8690857
commit b333dd0a96
8 changed files with 868 additions and 1 deletions

View File

@ -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

View File

@ -0,0 +1,9 @@
*.o
*.s
*.dump
*.aiff
*.aifc
*.table
/tabledesign_irix
/tabledesign_native
/tabledesign_orig

View File

@ -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

View File

@ -0,0 +1,104 @@
#include <stdlib.h>
#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);
}

View File

@ -0,0 +1,342 @@
#include <math.h>
#include <stdlib.h>
#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<j;i++) {
sum=a[i][j];
for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<=n;i++) {
sum=a[i][j];
for (k=1;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >= 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];
}
}

View File

@ -0,0 +1,89 @@
#include <stdio.h>
#include <stdlib.h>
#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;
}

View File

@ -0,0 +1,262 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <getopt.h>
#include <audiofile.h>
#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;
}

View File

@ -0,0 +1,30 @@
#ifndef TABLEDESIGN_H
#define TABLEDESIGN_H
#include <stdio.h>
#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