/**************************************************************************
 *                                                                        *
 * Copyright (c) 1996 Michael Richmond and Richard Treffers               *
 *                                                                        *
 *                    This software may be copied and distributed for     *
 *                    educational, research and not for profit services   *
 *                    provided that this copyright and statement are      *
 *                    included in all such copies.                        *
 *                                                                        *
 **************************************************************************/


/*	
	computes the Fourier transform of a series of number
		in the input file

	FFT1D filename 
		YCOL=		column number of 'y' data in file
		INVERSE		inverse transform
		BIAS		remove bias		(not default)
		NUM=		number of points to transform	
		STEP=		# of units between points (default=1.0)

	4/22/92 - now step is step -rrt
	9/22/92 - usage changed and call to exit(0) -rrt
			- get_col removed
			- stdout substituted for file and stdin allowed
	5/10/93 - help now calls exit
	4/22/94 - ournr.h instead of nr.h
	10/12/94 - NDATA increased
			 - add 'complex mode' output
	10/17/94 -skip lines containing # as first letter
	11/09/94 - num= now used 
	11/10/94 -  xcol dropped
				ycol defaults to 1
*/

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include "pcvista.h"
#include "ournr.h"

#define NDATA 262144	
#define NBUF	80
#ifndef PI
#define PI		3.14159265
#endif

#ifdef PROTO
void main(int, char **);
#else
#endif

#define USAGE \
  "usage: fft1d file [ycol=] [icol=] [bias] [inverse] [num=] [step=] \n      [complex] [verbose] [help]"

#define POLAR_MODE		0
#define COMPLEX_MODE	1

void 
main(argc, argv)
int argc;
char *argv[];
{
	int i, j, ndata, direction = 1;
	int  r_column = 1, number = 0, bias_flag = 0, gotfile = 0;
	int i_column=0;
	int mode=POLAR_MODE;
	int verbose = 0;
	FILE *fin;
	double mean, sum, sample_time = 1.0, freq_step, maxphase;
	double x, y, amp, freq, phase;
	static double data[2*NDATA], maxval, maxfreq;
	char buf[NBUF], fname[NBUF], *gotit;

	if(argc ==1){
		fprintf(stderr,"%s\n",USAGE);
		exit(0);
	}

	/* placate compiler */
	fin = NULL;

	for (i = 1; i < argc; i++) {
		if ((gotit = find("ycol", argv[i])) != NULL) {
			r_column = (int)(evaluate(gotit) + 0.5);
			continue;
		}
		if ((gotit = find("icol", argv[i])) != NULL) {
			i_column = (int)(evaluate(gotit) + 0.5);
			continue;
		}
		if ((gotit = find("num", argv[i]))  != NULL){
			number = (int)(evaluate(gotit) + 0.5);
			if (number > NDATA) {
				fprintf(stderr,"fft1d: maximum data points is %d\n", NDATA);
				exit(1);
			}
			continue;
		}
		if (keyword("inverse", argv[i])) {
			direction = 1;
			continue;
		}
		if (keyword("bias", argv[i])) {
			bias_flag = 1;
			continue;
		}
		if ((gotit = find("step", argv[i])) != NULL) {
			sample_time = evaluate(gotit);
			if (sample_time == 0.0)
				error(-1, "sample time cannot be zero");
			continue;
		}
		if (keyword("complex", argv[i])) {
			mode = COMPLEX_MODE;
			continue;
		}
		if (keyword("verbose", argv[i])) {
			verbose = 1;
			continue;
		}
		if (keyword("help", argv[i])) {
			fprintf(stderr,"%s\n",USAGE);
			exit(0);
		}
		if (gotfile) {
			fprintf(stderr,"fft1d: unknown option '%s'\n", argv[i]);
			exit(1);
		}
		strcpy(fname, argv[i]);
		gotfile = 1;
	}

	if (gotfile) {
		if(strcmp(fname,"-") == 0){
			fin = stdin;
		} else {
			if ((fin = fopen(fname, "r")) == NULL){
				fprintf(stderr, "fft1d:can't open input file '%s'\n",fname);
				exit(1);
			}	
		}
	}else{
		error(1,"fft1d: no file specified");
	}
	
	sum = 0.0;
	for (i = 0, j = 0; fgets(buf, NBUF, fin) != NULL; i++, j ++) {
		if(buf[0] == '#')	/* skip comment characters */
			continue;
		data[j] = get_col(buf,r_column);
		sum += data[j];
		j++;
		if(i_column)		/* fill in the imaginary data */
			data[j] = get_col(buf,i_column);
		else
			data[j]=0.;	
		if (i == NDATA) {
			fprintf(stderr, 
				"fft1d: number of data points can't exceed %d\n",NDATA);
			exit(1);
		}			
	}

	ndata = i;
	mean = sum/ndata;

	if (bias_flag) {
		if(mode == COMPLEX_MODE)
			fprintf(stderr,
				"fft1d: warning bias value removed from real part only\n");
		if (verbose)
			fprintf(stderr,"removing %8.3f from real data\n", mean);
		for (i = 0; i < ndata; i++)
			data[2*i] -= mean;		
	}

	if (number == 0) {
		for (number = 1; number<ndata; number *= 2)	/* find next power of two */
			;
	}
	else {
		if (ndata > number)							/* zero extra data */
			for (i = number; i < ndata; i++)
				data[2*i] = 0.0;		
	}
	if (ndata == 0)
		error(-1, "fft1d: no data points");
	if (verbose)
		fprintf(stderr,"ndata=%d  number=%d\n", ndata, number);

	four1(data - 1, number, direction);			/* this routine [1..2*n] */

	freq_step = 1.0/sample_time/number;

	if(verbose){
		if(mode == POLAR_MODE)
			printf("     Frequency    Amp.     Phase\n");
		if(mode == COMPLEX_MODE)
			printf("     Frequency    Real     Imaginary\n");
	}

	for (i = 0; i <= number; i += 2) {
		freq = i*freq_step/2.0;
		x = data[i];
		y = data[i + 1];
		if (i == number)
			x = x/2;
		switch(mode){
			case POLAR_MODE:
				amp = sqrt(x*x + y*y);
				if ((y != 0.0) || (x != 0.0))
					phase = 180*atan2(y, x)/PI;
				else
					phase = 0.0;
				if (amp > maxval) {
					maxval = amp;
					maxfreq = freq;
					maxphase = phase;
				}
				printf("%12.3f %10.6f %10.2f\n", freq, amp, phase);
				break;
			case COMPLEX_MODE:
				printf("%12.3f %10.6f %10.6f\n", freq, x, y);
				break;
			default:
				error(1,"fft1d: program error illegal mode");
		}
			
	}
	if (verbose && (mode=POLAR_MODE))
		fprintf(stderr,
		"Maximum abs val=%.3f at freq=%.4f and phase %.2f degs.\n",
			maxval, maxfreq, maxphase);
	exit(0);
}