soundmodem/fsk/genpsptbl.c

/*****************************************************************************/

/*
 *      genpsptbl.c  --  Per Survivor Processing Table Generator.
 *
 *      Copyright (C) 2000
 *        Thomas Sailer (sailer@ife.ee.ethz.ch)
 *
 *      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., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */

/*****************************************************************************/

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>
#include <stdio.h>

#include "mat.h"

/* --------------------------------------------------------------------- */

#define DEBUG

#define VECLENGTH  3

#define ESTDATASYMS 8
#define ESTPARAMS   (VECLENGTH+1)
#define ESTSYMBOLS  (ESTDATASYMS-VECLENGTH+1)
#define ESTSHIFT    8

/* ---------------------------------------------------------------------- */

static int fcmatprintf(FILE *f, const char *name, unsigned int size1, unsigned int stride1,
                       unsigned int size2, unsigned int stride2, const cplxfloat_t *m)
{
        unsigned int i, j;
        int ret = 0;

        fprintf(f, "%s = [", name);
        for (i = 0; i < size1; i++) {
                for (j = 0; j < size2; j++) {
                        ret += fprintf(f, " %g", real(m[i*stride1 + j*stride2]));
                        if (imag(m[i*stride1 + j*stride2]) != 0)
                                ret += fprintf(f, "%+gi", imag(m[i*stride1 + j*stride2]));
                }
                if (i+1 < size1)
                        ret += fprintf(f, " ; ...\n    ");
        }
        ret += fprintf(f, " ];\n");
        return ret;
}

/* --------------------------------------------------------------------- */

static void gentrainmats(FILE *f)
{
        cplxfloat_t ma[ESTSYMBOLS*ESTPARAMS], mah[ESTPARAMS*ESTSYMBOLS], maha[ESTPARAMS*ESTPARAMS];
        cplxfloat_t mahainv[ESTPARAMS*ESTPARAMS], mahainvah[ESTPARAMS*ESTSYMBOLS];
	cplxfloat_t det;
	unsigned int i, j, k;
	unsigned long singular[1<<(ESTDATASYMS-5)];

	fprintf(f, "/* this file is automatically generated, do not edit!! */\n\n"
		"#define VECLENGTH   %u\n"
		"#define ESTDATAMASK %u\n"
		"#define ESTDATASYMS %u\n"
		"#define ESTPARAMS   %u\n"
		"#define ESTSYMBOLS  %u\n"
		"#define ESTSHIFT    %u\n"
		"\n\nstatic const int16_t estmat[%u][ESTPARAMS*ESTSYMBOLS] = {",
		VECLENGTH, (1<<ESTDATASYMS)-1, ESTDATASYMS, ESTPARAMS, ESTSYMBOLS,
		ESTSHIFT, 1<<ESTDATASYMS);
	memset(singular, 0, sizeof(singular));
	for (i = 0; i < (1<<ESTDATASYMS); i++) {
		if (i)
			fputc(',', f);
		fprintf(f, "\n\t");
		for (j = 0; j < ESTSYMBOLS; j++) {
			cplx(ma[j*ESTPARAMS], 1, 0);
			for (k = 0; k < VECLENGTH; k++)
				cplx(ma[j*ESTPARAMS+k+1], (i & (1 << (j+VECLENGTH-1-k))) ? 1 : -1, 0);
		}
#ifdef DEBUG
		fprintf(f, "/*\n");
		fcmatprintf(f, "a", ESTSYMBOLS, ESTPARAMS, ESTPARAMS, 1, ma);
#endif
		/* transpose it */
		fchermtranspose(mah, ma, ESTSYMBOLS, ESTPARAMS);
		fcmul(maha, mah, ma, ESTPARAMS, ESTSYMBOLS, ESTPARAMS);
		det = fcdet(maha, ESTPARAMS);
		if (det.re == 0) {
			singular[i >> 5] |= 1 << (i & 31);
			for (j = 0; j < ESTPARAMS*ESTSYMBOLS; j++)
				cplx(mahainvah[j], 0, 0);
#ifdef DEBUG
			fprintf(f, " singular\n");
#endif
		} else {
			fcinv(mahainv, maha, ESTPARAMS);
			fcmul(mahainvah, mahainv, mah, ESTPARAMS, ESTPARAMS, ESTSYMBOLS);
#ifdef DEBUG
			fcmatprintf(f, "ahainvah", ESTPARAMS, ESTSYMBOLS, ESTSYMBOLS, 1, mahainvah);
#endif
		}
#ifdef DEBUG
		fprintf(f, "*/\n\t");
#endif
		fprintf(f, "{ ");
		for (j = 0;; j++) {
			fprintf(f, "%d", (int)(mahainvah[j].re * (double)(1<<ESTSHIFT)));
			if (j >= (ESTPARAMS*ESTSYMBOLS-1))
				break;
			fprintf(f, ", ");
		}
		fprintf(f, " }");
	}
	fprintf(f, "\n};\n\nstatic const unsigned long estsingular[%u] = {\n\t",
		1<<(ESTDATASYMS-5));
	for (i = 0; ; i++) {
		fprintf(f, "0x%08lx", singular[i]);
		if (i >= (1<<(ESTDATASYMS-5))-1)
			break;
		fprintf(f, ", ");
	}
	fprintf(f, "\n};\n\n");
}

/* ---------------------------------------------------------------------- */

int main(int argc, char *argv[])
{
	gentrainmats(stdout);
        return 0;
}