#ifndef lint
static char sccsid[] = "@(#)fftanal.c	1.4 (mg@ukc.ac.uk) 6/15/90";
#endif
/*
 *	Show power spectrum of a signal over time as a greyscaled rectangle.
 *	X-axis is time, y-axis is frequency, darkness represents power.
 *
 * Input is a file of single-precision floats representing real point data,
 * Output is a sun raster file.
 *
 * Usage: fftanal [-wwindow_size] [inputfilename]
 */

#include <stdio.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <math.h>	/* for atof */
#include "fft.h"

extern char *mmalloc();		/* csound's malloc that checks for failure */

#define MAX_COLOUR 31	/* max value for grey-level quantisation of image */

static int srate = 0;
static int window_type = HANNING;
static int window_size = 0;
static float period = 0.0;
static char *infilename = NULL;
static FILE *infp = NULL;

static char *progname;

/* Stats gathering */
static float max_max = 0.0;

main(argc, argv)
char **argv;
{
	int i;

	progname = argv[0];

	/* Process arguments */
	for (i=1; i<argc;i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
			case 'w':	/* Window size.
					 * (twice the no of output points) */
				window_size = atoi(&argv[i][2]);
				break;
			case 's':	/* sampling rate */
				srate = atoi(&argv[i][2]);
				break;
			case 'p':	/* period with which we do ffts */
				period = atof(&argv[i][2]);
				break;
			default:
				fprintf(stderr, "%s: Unknown option \"%s\".\n",
					progname, argv[i]);
				exit(1);
			}
		} else {
			infp = fopen(argv[i], "r");
			if (infp == NULL) {
				fprintf(stderr, "%s: ", argv[0]);
				perror(argv[i]);
				exit(1);
			}
		}
	}

	if (infp == NULL) {
		fputs("Input file name is compulsory\n", stderr);
		exit(1);
	}

	/* Supply defaults */
	if (srate == 0) srate = 8192;
	if (window_size == 0) window_size = 1024;
	if (period == 0.0) period = ((double)window_size)/srate;

	process(infp, stdout);

	exit(0);
}

process(infp, outfp)
FILE *infp;	/* handle for file of floats */
FILE *outfp;	/* handle for sun raster output file */
{
	float *inbuf;		/* buffer of input values */
	float *outbuf;		/* buffer of output values */
	int ncoefs = window_size / 2;	/* size of input buffer */
	int nsamples;		/* Number of samples in the file */
	int nframes;			/* number of fft frames */
	int i;

	/* Find size of file */
	{
		struct stat stbuf;

		if (fstat(fileno(infp), &stbuf) != 0) {
			fprintf(stderr, "%s: fstat failed on ", progname);
			perror(infilename);
			exit(1);
		}

		nsamples = stbuf.st_size / sizeof(float);
	}

	/* see how many frames this will give us at the given rates */
	if (window_size > nsamples) {
		fprintf(stderr, "%s: Not enough samples for a single fft.\n",
			progname);
		exit(1);
	}
	nframes = (nsamples - window_size)/(period * srate) + 1;

	fprintf(stderr, "nsamples = %ld, window_size = %ld, period = %g, srate = %ld, nframes = %ld\n",
		nsamples, window_size, period, srate, nframes);

	/* Get collection buffer for sample data */
	inbuf = (float *) mmalloc((long) (ncoefs * 2) * sizeof(float));
	/* and for FFT results */
	outbuf = (float *) mmalloc((long) ncoefs * sizeof(float));

	sunrf_open(ncoefs, nframes, (unsigned char)MAX_COLOUR, outfp);

	for (i=0; i<nframes; i++) {
		get_samples(period*i, inbuf, window_size);

		/* process one frame */
		fft(FORWARD, window_size, window_type,
		    inbuf, REAL, LINEAR,
		    outbuf, MAG, LINEAR);
		
		display_line(outbuf, ncoefs, outfp);
	}

	sunrf_close(outfp);
}

/* Read N samples from input file into buffer */
get_samples(start_time, buf, n)
float start_time;	/* time at which to start reading */
float *buf;
{
	int nread;
	long offset;

	/* find where to start reading, ensuring alignment. */
	/* "When a floating type is converted to an integral type, the
	 * fractional part is discarded", (K&R 2) */
	offset = ((long)(start_time * srate)) * sizeof(float);

	if (fseek(infp, offset, 0) != 0) {
		fprintf(stderr, "%s: Cannot seek to position %ld on input file.\n",
			progname, offset);
		exit(1);
	}

	nread = fread((char *)buf, sizeof(float), n, infp);
	if (nread == 0 && ferror(infp)) {
		fprintf(stderr, "%s: Read error on ", progname);
		perror(infilename);
		exit(1);
	}
		
	return(nread);
}

display_line(buf, ncoefs, ofp)
float *buf;	/* raw data */
int ncoefs;	/* line width */
FILE *ofp;	/* output file pointer */
{
	float max;
	float *fp;
	static unsigned char *pbuf = NULL;
	unsigned char *pp;
	register int i;
	register float sample;
	int errcnt = 0;

	/* convert amplitude to power and find maximum pixel value */
	for (i=0, fp=buf, max=0.0; i<ncoefs; i++, fp++) {
		sample = (*fp = *fp * i);
		if (sample > max) max = sample;
		if (sample < 0) {
			errcnt++;
			*fp = max;
		}
	}

	if (errcnt) {
		fprintf(stderr, "%s: %ld negative output values\n",
			progname, errcnt);
	}

	/* avoid bottoming out */
	if (max == 0.0) {
		max = 1.0;
	}

	if (max > max_max) max_max = max;

	/* ensure pixel buffer is allocated */
	if (pbuf == NULL) pbuf = (unsigned char *) mmalloc((long) ncoefs);

	/* convert from float values to pixel values */
	for (i=0, fp=buf, pp=pbuf; i<ncoefs; i++) {
		*pp++ = *fp++ * MAX_COLOUR / max_max;
	}

	/* and output it */
	sunrf_line(pbuf, ncoefs, stdout);
}