/* dft.c -- calculates discrete Fourier transform of binary short sound samples
/*  (from stdin) and writes binary float coefficients (mag or db) to stdout.
/* usage:	dft [-hdp] windowsiz stepsiz  <infile  >outfile
/*		     -h for hanning window, -d for dbout, -p for printout
/*		defaults are rectangular window, magnitude coeffs, noprint
/*		windowsiz: suggest equiv of 5 pitch periods, in samples
/*		stepsiz:   a bit less for overlapped frames		*/

#include "sysdep.h"
#include <stdio.h>
#include <math.h>

#define	WINDMAX	500

short	sampbuf[WINDMAX];	   	/* buffer for input samples	   */
float	outbuf[WINDMAX/2];		/* frame of out coefs in mag or db */
double	hanwindow[WINDMAX];		/* hanning window values	   */

main(argc,argv)
 int argc;
 char **argv;
{
	double	theta, a, b, c, db, *wsin, *wcos;
	double	pidws, twopidws, twodws, n2pidws;
	int	windsiz, stepsiz, windbytes, stepbytes, ncoefs;
	int	hanning = 0, dbout = 0, print = 0, frmcnt = 0;
	char	*cp;
register int	n, k;
register short	*samp, *samp2;
register double *sinp, *cosp, *hanp;

	argc--;	argv++;				/* align onto command args */
	while (argc && (cp = *argv) && *cp++ == '-') {
		while (*cp != '\0')     	/* for all flag arguments: */
			switch(*cp++) {
			case 'h': hanning = 1;	/*	-h, hanning window */
				  break;
			case 'd': dbout = 1;	/*	-d, output is db */
				  break;
			case 'p': print = 1;	/*	-p, set printout */
				  break;
			default:  die("unknown flag");
			}
		argc--;	argv++;
	}
	if (argc < 2)
		die("usage: dft [-hdp] windsiz stepsiz <infile >outfile");
	sscanf(*argv++, "%d", &windsiz);	/* rd window size, step size */
	sscanf(*argv++, "%d", &stepsiz);
	if (windsiz > WINDMAX)
		die("windsiz too large for program;  recompile");
	windbytes = windsiz * 2;
	stepbytes = stepsiz * 2;
	ncoefs = windsiz / 2;
	twodws = 2.0 / windsiz;
	twopidws = 6.2831854 / windsiz;
	pidws = 3.1415927 / windsiz;
	if (hanning) {				/* if hanning requested,    */
		hanp = hanwindow;
		for (k=1; k<=windsiz; k++) {	/*    for this window size  */
			a = sin(k*pidws);	/*    create hanning points */
			*hanp++ = a * a;
		}
	}
	fprintf(stderr,"dft: %s window, %s out, making tables ...\n",
		(hanning) ? "hanning":"rect", (dbout) ? "db":"magnitude");
	wsin = (double *)mmalloc((long)ncoefs*windsiz*8); /* alloc two 2-dim */
	wcos = (double *)mmalloc((long)ncoefs*windsiz*8); /* weightd sintabs */
	for (n=1,sinp=wsin,cosp=wcos; n<=ncoefs; n++) {   /* now fill tables */
		n2pidws = n * twopidws;
		if (hanning) {
			hanp = hanwindow;
			for (k=1; k<=windsiz; k++) {	  /*   with sines    */
				theta = k * n2pidws;
				a = twodws * *hanp++;	  /*   times hanning */
				*sinp++ = a * sin(theta);
				*cosp++ = a * cos(theta);
			}
		}
		else {
			for (k=1; k<=windsiz; k++) {	   /* or sines alone */
				theta = k * n2pidws;	   /*  (rect window) */
				*sinp++ = twodws*sin(theta);
				*cosp++ = twodws*cos(theta);
			}
		}
	}
	n = read(0, (char *)sampbuf, windbytes);  /* init samp window w. input */
	if (n != windbytes)
		die("premature end of infile");

frame:	for (n=1,sinp=wsin,cosp=wcos; n<=ncoefs; n++) {	   /* for one frame: */
		a = 0.0;
		b = 0.0;
		samp = sampbuf;
		for (k=1; k<=windsiz; k++) {		/*  calculate coefs  */
			a += *samp * *cosp++;
			b += *samp++ * *sinp++;
		}
		c = sqrt( a*a + b*b );
		if (!dbout)
			outbuf[n-1] = c;		/* stor as magnitude */
		else outbuf[n-1] = db = 20. * log10(c);	/*	or db	     */
		if (print)
			fprintf(stderr,"num %d coef %f db %f per %d freq %d\n",
				n, c, db, windsiz/n, 32000*n/windsiz );
	}
	write(1, (char *)outbuf, ncoefs*4);	/*  then wrt frame to stdout */
	fprintf(stderr,"dft frame %d\n",++frmcnt);

	if ((n = windsiz - stepsiz) > 0) {	/* nxt: for step < windowsiz */
		samp = sampbuf;
		samp2 = sampbuf + stepsiz;
		while (n--)
			*samp++ = *samp2++;	/*	slide the sample buf */
		n = read(0, (char *)samp, stepbytes); /* & refill to windsiz */
		if (n == stepbytes)
			goto frame;
	}
	else {
		for (n = -n; n>=windsiz; n-=windsiz)	/* else waste any    */
			read(0, (char *)sampbuf, windbytes); /*  extra samps */
		if (n)	read(0, (char *)sampbuf, n*2);
		n = read(0,(char *)sampbuf, windbytes); /* then rd whole new window */
		if (n == windbytes);
			goto frame;		        /*	& go calc new frame */
	}
	fprintf(stderr,"wrote %d frames of %d coefs\n", frmcnt, ncoefs);
}


die(s)
 char *s;
{
	fprintf(stderr,"dft: %s\n",s);
	exit(1);
}