#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifndef __STDC__
#  ifndef const
#    define const
#  endif
#endif

/* This tells Alpha OSF/1 not to define a getopt prototype in <stdio.h>.  */
#ifndef _NO_PROTO
#define _NO_PROTO
#endif

#include <stdio.h>
#include "tailor.h"

#if defined (_LIBC) || !defined (__GNU_LIBRARY__)

#ifdef	__GNU_LIBRARY__

#include <stdlib.h>
#endif	/* GNU C library.  */

#include "getopt.h"

char *optarg = 0;

/* XXX 1003.2 says this must be 1 before any call.  */
int optind = 0;

static char *nextchar;

int opterr = 1;

#define BAD_OPTION '\0'
int optopt = BAD_OPTION;

static enum
{
  REQUIRE_ORDER, PERMUTE, RETURN_IN_ORDER
} ordering;

#ifdef	__GNU_LIBRARY__

#include <string.h>
#define	my_index	strchr
#define	my_strlen	strlen
#else

/* Avoid depending on library functions or files
   whose names are inconsistent.  */

#if __STDC__ || defined(PROTO)
extern char *getenv(const char *name);
extern int  strcmp (const char *s1, const char *s2);
extern int  strncmp(const char *s1, const char *s2, int n);

static int my_strlen(const char *s);
static char *my_index (const char *str, int chr);
#else
extern char *getenv ();
#endif

static int
my_strlen (str)
     const char *str;
{
  int n = 0;
  while (*str++)
    n++;
  return n;
}

static char *
my_index (str, chr)
     const char *str;
     int chr;
{
  while (*str)
    {
      if (*str == chr)
	return (char *) str;
      str++;
    }
  return 0;
}

#endif				/* GNU C library.  */

/* Handle permutation of arguments.  */

static int first_nonopt;
static int last_nonopt;

#if __STDC__ || defined(PROTO)
static void exchange (char **argv);
#endif

static void
exchange (argv)
     char **argv;
{
  char *temp, **first, **last;

  /* Reverse all the elements [first_nonopt, optind) */
  first = &argv[first_nonopt];
  last  = &argv[optind-1];
  while (first < last) {
    temp = *first; *first = *last; *last = temp; first++; last--;
  }
  /* Put back the options in order */
  first = &argv[first_nonopt];
  first_nonopt += (optind - last_nonopt);
  last  = &argv[first_nonopt - 1];
  while (first < last) {
    temp = *first; *first = *last; *last = temp; first++; last--;
  }

  /* Put back the non options in order */
  first = &argv[first_nonopt];
  last_nonopt = optind;
  last  = &argv[last_nonopt-1];
  while (first < last) {
    temp = *first; *first = *last; *last = temp; first++; last--;
  }
}

int
_getopt_internal (argc, argv, optstring, longopts, longind, long_only)
     int argc;
     char *const *argv;
     const char *optstring;
     const struct option *longopts;
     int *longind;
     int long_only;
{
  int option_index;

  optarg = 0;

  if (optind == 0)
    {
      first_nonopt = last_nonopt = optind = 1;

      nextchar = NULL;

      /* Determine how to handle the ordering of options and nonoptions.  */

      if (optstring[0] == '-')
	{
	  ordering = RETURN_IN_ORDER;
	  ++optstring;
	}
      else if (optstring[0] == '+')
	{
	  ordering = REQUIRE_ORDER;
	  ++optstring;
	}
      else if (getenv ("POSIXLY_CORRECT") != NULL)
	ordering = REQUIRE_ORDER;
      else
	ordering = PERMUTE;
    }

  if (nextchar == NULL || *nextchar == '\0')
    {
      if (ordering == PERMUTE)
	{
	  /* If we have just processed some options following some non-options,
	     exchange them so that the options come first.  */

	  if (first_nonopt != last_nonopt && last_nonopt != optind)
	    exchange ((char **) argv);
	  else if (last_nonopt != optind)
	    first_nonopt = optind;

	  /* Now skip any additional non-options
	     and extend the range of non-options previously skipped.  */

	  while (optind < argc
		 && (argv[optind][0] != '-' || argv[optind][1] == '\0')
#ifdef GETOPT_COMPAT
		 && (longopts == NULL
		     || argv[optind][0] != '+' || argv[optind][1] == '\0')
#endif				/* GETOPT_COMPAT */
		 )
	    optind++;
	  last_nonopt = optind;
	}

          if (optind != argc && !strcmp (argv[optind], "--"))
	{
	  optind++;

	  if (first_nonopt != last_nonopt && last_nonopt != optind)
	    exchange ((char **) argv);
	  else if (first_nonopt == last_nonopt)
	    first_nonopt = optind;
	  last_nonopt = argc;

	  optind = argc;
	}

   

      if (optind == argc)
	{
	  /* Set the next-arg-index to point at the non-options
	     that we previously skipped, so the caller will digest them.  */
	  if (first_nonopt != last_nonopt)
	    optind = first_nonopt;
	  return EOF;
	}

      /* If we have come to a non-option and did not permute it,
	 either stop the scan or describe it to the caller and pass it by.  */

      if ((argv[optind][0] != '-' || argv[optind][1] == '\0')
#ifdef GETOPT_COMPAT
	  && (longopts == NULL
	      || argv[optind][0] != '+' || argv[optind][1] == '\0')
#endif				/* GETOPT_COMPAT */
	  )
	{
	  if (ordering == REQUIRE_ORDER)
	    return EOF;
	  optarg = argv[optind++];
	  return 1;
	}

      /* We have found another option-ARGV-element.
	 Start decoding its characters.  */

      nextchar = (argv[optind] + 1
		  + (longopts != NULL && argv[optind][1] == '-'));
    }

  if (longopts != NULL
      && ((argv[optind][0] == '-'
	   && (argv[optind][1] == '-' || long_only))
#ifdef GETOPT_COMPAT
	  || argv[optind][0] == '+'
#endif				/* GETOPT_COMPAT */
	  ))
    {
      const struct option *p;
      char *s = nextchar;
      int exact = 0;
      int ambig = 0;
      const struct option *pfound = NULL;
      int indfound = 0;

      while (*s && *s != '=')
	s++;

      /* Test all options for either exact match or abbreviated matches.  */
      for (p = longopts, option_index = 0; p->name;
	   p++, option_index++)
	if (!strncmp (p->name, nextchar, s - nextchar))
	  {
	    if (s - nextchar == my_strlen (p->name))
	      {
		/* Exact match found.  */
		pfound = p;
		indfound = option_index;
		exact = 1;
		break;
	      }
	    else if (pfound == NULL)
	      {
		/* First nonexact match found.  */
		pfound = p;
		indfound = option_index;
	      }
	    else
	      /* Second nonexact match found.  */
	      ambig = 1;
	  }

      if (ambig && !exact)
	{
	  if (opterr)
	    fprintf (stderr, "%s: option `%s' is ambiguous\n",
		     argv[0], argv[optind]);
	  nextchar += my_strlen (nextchar);
	  optind++;
	  return BAD_OPTION;
	}

      if (pfound != NULL)
	{
	  option_index = indfound;
	  optind++;
	  if (*s)
	    {
	      /* Don't test has_arg with >, because some C compilers don't
		 allow it to be used on enums.  */
	      if (pfound->has_arg)
		optarg = s + 1;
	      else
		{
		  if (opterr)
		    {
		      if (argv[optind - 1][1] == '-')
			/* --option */
			fprintf (stderr,
				 "%s: option `--%s' doesn't allow an argument\n",
				 argv[0], pfound->name);
		      else
			/* +option or -option */
			fprintf (stderr,
			     "%s: option `%c%s' doesn't allow an argument\n",
			     argv[0], argv[optind - 1][0], pfound->name);
		    }
		  nextchar += my_strlen (nextchar);
		  return BAD_OPTION;
		}
	    }
	  else if (pfound->has_arg == 1)
	    {
	      if (optind < argc)
		optarg = argv[optind++];
	      else
		{
		  if (opterr)
		    fprintf (stderr, "%s: option `%s' requires an argument\n",
			     argv[0], argv[optind - 1]);
		  nextchar += my_strlen (nextchar);
		  return optstring[0] == ':' ? ':' : BAD_OPTION;
		}
	    }
	  nextchar += my_strlen (nextchar);
	  if (longind != NULL)
	    *longind = option_index;
	  if (pfound->flag)
	    {
	      *(pfound->flag) = pfound->val;
	      return 0;
	    }
	  return pfound->val;
	}

      if (!long_only || argv[optind][1] == '-'
#ifdef GETOPT_COMPAT
	  || argv[optind][0] == '+'
#endif				/* GETOPT_COMPAT */
	  || my_index (optstring, *nextchar) == NULL)
	{
	  if (opterr)
	    {
	      if (argv[optind][1] == '-')
		/* --option */
		fprintf (stderr, "%s: unrecognized option `--%s'\n",
			 argv[0], nextchar);
	      else
		/* +option or -option */
		fprintf (stderr, "%s: unrecognized option `%c%s'\n",
			 argv[0], argv[optind][0], nextchar);
	    }
	  nextchar = (char *) "";
	  optind++;
	  return BAD_OPTION;
	}
    }

  /* Look at and handle the next option-character.  */

  {
    char c = *nextchar++;
    char *temp = my_index (optstring, c);

    /* Increment `optind' when we start to process its last character.  */
    if (*nextchar == '\0')
      ++optind;

    if (temp == NULL || c == ':')
      {
	if (opterr)
	  {
#if 0
	    if (c < 040 || c >= 0177)
	      fprintf (stderr, "%s: unrecognized option, character code 0%o\n",
		       argv[0], c);
	    else
	      fprintf (stderr, "%s: unrecognized option `-%c'\n", argv[0], c);
#else
	    /* 1003.2 specifies the format of this message.  */
	    fprintf (stderr, "%s: illegal option -- %c\n", argv[0], c);
#endif
	  }
	optopt = c;
	return BAD_OPTION;
      }
    if (temp[1] == ':')
      {
	if (temp[2] == ':')
	  {
	    /* This is an option that accepts an argument optionally.  */
	    if (*nextchar != '\0')
	      {
		optarg = nextchar;
		optind++;
	      }
	    else
	      optarg = 0;
	    nextchar = NULL;
	  }
	else
	  {
	    /* This is an option that requires an argument.  */
	    if (*nextchar != '\0')
	      {
		optarg = nextchar;
		/* If we end this ARGV-element by taking the rest as an arg,
		   we must advance to the next element now.  */
		optind++;
	      }
	    else if (optind == argc)
	      {
		if (opterr)
		  {
#if 0
		    fprintf (stderr, "%s: option `-%c' requires an argument\n",
			     argv[0], c);
#else
		    /* 1003.2 specifies the format of this message.  */
		    fprintf (stderr, "%s: option requires an argument -- %c\n",
			     argv[0], c);
#endif
		  }
		optopt = c;
		if (optstring[0] == ':')
		  c = ':';
		else
		  c = BAD_OPTION;
	      }
	    else
	      /* We already incremented `optind' once;
		 increment it again when taking next ARGV-elt as argument.  */
	      optarg = argv[optind++];
	    nextchar = NULL;
	  }
      }
    return c;
  }
}

int
getopt (argc, argv, optstring)
     int argc;
     char *const *argv;
     const char *optstring;
{
  return _getopt_internal (argc, argv, optstring,
			   (const struct option *) 0,
			   (int *) 0,
			   0);
}

int
getopt_long (argc, argv, options, long_options, opt_index)
     int argc;
     char *const *argv;
     const char *options;
     const struct option *long_options;
     int *opt_index;
{
  return _getopt_internal (argc, argv, options, long_options, opt_index, 0);
}

#endif	/* _LIBC or not __GNU_LIBRARY__.  */

#ifdef TEST

/* Compile with -DTEST to make an executable for use in testing
   the above definition of `getopt'.  */

int
main (argc, argv)
     int argc;
     char **argv;
{
  int c;
  int digit_optind = 0;

  while (1)
    {
      int this_option_optind = optind ? optind : 1;

      c = getopt (argc, argv, "abc:d:0123456789");
      if (c == EOF)
	break;

      switch (c)
	{
	case '0':
	case '1':
	case '2':
	case '3':
	case '4':
	case '5':
	case '6':
	case '7':
	case '8':
	case '9':
	  if (digit_optind != 0 && digit_optind != this_option_optind)
	    printf ("digits occur in two different argv-elements.\n");
	  digit_optind = this_option_optind;
	  printf ("option %c\n", c);
	  break;

	case 'a':
	  printf ("option a\n");
	  break;

	case 'b':
	  printf ("option b\n");
	  break;

	case 'c':
	  printf ("option c with value `%s'\n", optarg);
	  break;

	case BAD_OPTION:
	  break;

	default:
	  printf ("?? getopt returned character code 0%o ??\n", c);
	}
    }

  if (optind < argc)
    {
      printf ("non-option ARGV-elements: ");
      while (optind < argc)
	printf ("%s ", argv[optind++]);
      printf ("\n");
    }

  exit (0);
}

#endif /* TEST */

#include <stdio.h>
#include <math.h>
#ifdef __STDC__
#include <stdlib.h>
#else
extern void exit();
#endif

#ifndef BSD
# include <string.h>
#else           /* for Berkeley UNIX only */
# include <strings.h>
# define strchr index
#endif

#define	LEN	16384		/* maximum points in FFT */
#ifdef i386
#define strcasecmp strcmp
#endif
#define PI	M_PI	/* pi to machine precision, defined in math.h */
#define TWOPI	(2.0*PI)

void four1();
void realft();

double wsum;

/* See Oppenheim & Schafer, Digital Signal Processing, p. 241 (1st ed.) */
double win_bartlett(j, n)
int j, n;
{
    double a = 2.0/(n-1), w;

    if ((w = j*a) > 1.0) w = 2.0 - w;
    wsum += w;
    return (w);
}

/* See Oppenheim & Schafer, Digital Signal Processing, p. 242 (1st ed.) */
double win_blackman(j, n)
int j, n;
{
    double a = 2.0*PI/(n-1), w;

    w = 0.42 - 0.5*cos(a*j) + 0.08*cos(2*a*j);
    wsum += w;
    return (w);
}

/* See Harris, F.J., "On the use of windows for harmonic analysis with the
   discrete Fourier transform", Proc. IEEE, Jan. 1978 */
double win_blackman_harris(j, n)
int j, n;
{
    double a = 2.0*PI/(n-1), w;

    w = 0.35875 - 0.48829*cos(a*j) + 0.14128*cos(2*a*j) - 0.01168*cos(3*a*j);
    wsum += w;
    return (w);
}

/* See Oppenheim & Schafer, Digital Signal Processing, p. 242 (1st ed.) */
double win_hamming(j, n)
int j, n;
{
    double a = 2.0*PI/(n-1), w;

    w = 0.54 - 0.46*cos(a*j);
    wsum += w;
    return (w);
}

/* See Oppenheim & Schafer, Digital Signal Processing, p. 242 (1st ed.)
   The second edition of Numerical Recipes calls this the "Hann" window. */
double win_hanning(j, n)
int j, n;
{
    double a = 2.0*PI/(n-1), w;

    w = 0.5 - 0.5*cos(a*j);
    wsum += w;
    return (w);
}

/* See Press, Flannery, Teukolsky, & Vetterling, Numerical Recipes in C,
   p. 442 (1st ed.) */
double win_parzen(j, n)
int j, n;
{
    double a = (n-1)/2.0, w;

    if ((w = (j-a)/(a+1)) > 0.0) w = 1 - w;
    else w = 1 + w;
    wsum += w;
    return (w);
}

/* See any of the above references. */
double win_square(j, n)
int j, n;
{
    wsum += 1.0;
    return (1.0);
}

/* See Press, Flannery, Teukolsky, & Vetterling, Numerical Recipes in C,
   p. 442 (1st ed.) or p. 554 (2nd ed.) */
double win_welch(j, n)
int j, n;
{
    double a = (n-1)/2.0, w;

    w = (j-a)/(a+1);
    w = 1 - w*w;
    wsum += w;
    return (w);
}

char *pname;
double rsum;
FILE *ifile;
int m, n;
int cflag;
int decimation = 1;
int iflag;
int fflag;
int len = LEN;
int nflag;
int Nflag;
int nchunks;
int pflag;
int Pflag;
int smooth = 1;
int wflag;
int zflag;
static float *c;
double freq, fstep, norm, rmean, (*window)();

main(argc, argv)
int argc;
char *argv[];
{
    int i;
    char *prog_name();
    double atof();
    void help();

    pname = prog_name(argv[0]);
    if (--argc < 1) {
	help();
	exit(1);
    }
    else if (strcmp(argv[argc], "-") == 0)
	ifile = stdin;		/* read data from standard input */
    else if ((ifile = fopen(argv[argc], "rt")) == NULL) {
	fprintf(stderr, "%s: can't open %s\n", pname, argv[argc]);
	exit(2);
    }
    for (i = 1; i < argc; i++) {
	if (*argv[i] == '-') switch (*(argv[i]+1)) {
	  case 'c':	/* print complex FFT (rectangular form) */
	    cflag = 1;
	    break;
	  case 'f':	/* print frequencies */
	    if (++i >= argc) {
		fprintf(stderr, "%s: sampling frequency must follow -f\n",
			pname);
		exit(1);
	    }
	    freq = atof(argv[i]);
	    fflag = 1;
	    break;
	  case 'h':	/* print help and exit */
	    help();
	    exit(0);
	    break;
	  case 'i':	/* calculate inverse FFT from `fft -c' format input */
	    if (argc > 2) {
		fprintf(stderr, "%s: no other option may be used with -i\n",
			pname);
		exit(1);
	    }
	    iflag = 1;
	    break;
	  case 'I':	/* calculate inverse FFT from `fft -p' format input */
	    if (argc > 2) {
		fprintf(stderr, "%s: no other option may be used with -I\n",
			pname);
		exit(1);
	    }
	    iflag = -1;
	    break;
	  case 'l':	/* perform up to n-point transforms */
	    if (++i >= argc) {
		fprintf(stderr, "%s: transform size must follow -l\n", pname);
		exit(1);
	    }
	    len = atoi(argv[i]);
	    break;
	  case 'n':	/* process in overlapping n-point chunks, output avg */
	    if (++i >= argc) {
		fprintf(stderr, "%s: chunk size must follow -n\n", pname);
		exit(1);
	    }
	    nflag = atoi(argv[i]);
	    break;
	  case 'N':	/* process in overlapping n-point chunks, output raw */
	    if (++i >= argc) {
		fprintf(stderr, "%s: chunk size must follow -N\n", pname);
		exit(1);
	    }
	    Nflag = atoi(argv[i]);
	    break;
	  case 'p':	/* print phases (polar form) */
	    pflag = 1;
	    break;
	  case 'P':	/* print power spectrum (squared magnitudes) */
	    Pflag = 1;
	    break;
	  case 's':	/* smooth spectrum */
	    if (++i >= argc || ((smooth=atoi(argv[i])) < 2) || smooth > 1024) {
		fprintf(stderr, "%s: smoothing parameter must follow -s\n",
			pname);
		exit(1);
	    }
	    break;
	  case 'S':	/* smooth and decimate spectrum */
	    if (++i >= argc || ((smooth=atoi(argv[i])) < 2) || smooth > 1024) {
		fprintf(stderr, "%s: smoothing parameter must follow -s\n",
			pname);
		exit(1);
	    }
	    decimation = smooth;
	    break;
	  case 'w':	/* apply windowing function to input */
	    if (++i >= argc) {
		fprintf(stderr, "%s: window type must follow -w\n",
			pname);
		exit(1);
	    }
	    if (strcasecmp(argv[i], "Bartlett") == 0)
		window = win_bartlett;
	    else if (strcasecmp(argv[i], "Blackman") == 0)
		window = win_blackman;
	    else if (strcasecmp(argv[i], "Blackman-Harris") == 0)
		window = win_blackman_harris;
	    else if (strcasecmp(argv[i], "Hamming") == 0)
		window = win_hamming;
	    /* Numerical Recipes 2nd ed. calls Hanning window "Hann window" */
	    else if (strcasecmp(argv[i], "Hann") == 0)
		window = win_hanning;
	    else if (strcasecmp(argv[i], "Hanning") == 0)
		window = win_hanning;
	    else if (strcasecmp(argv[i], "Parzen") == 0)
		window = win_parzen;
	    else if (strcasecmp(argv[i], "Square") == 0 ||
		     strcasecmp(argv[i], "Rectangular") == 0 ||
		     strcasecmp(argv[i], "Dirichlet") == 0)
		window = win_square;
	    else if (strcasecmp(argv[i], "Welch") == 0)
		window = win_welch;
	    else {
		fprintf(stderr, "%s: unrecognized window type %s\n",
			pname, argv[i]);
		exit(1);
	    }
	    wflag = 1;
	    break;
	  case 'z':	/* zero-mean the input */
	    zflag = 1;
	    break;
	  case 'Z':	/* zero-mean and detrend the input */
	    zflag = 2;
	    break;
	  default:
	    fprintf(stderr, "%s: unrecognized option %s ignored\n",
		    pname, argv[i]);
	    break;
	}
    }
    if (cflag) {
	if (fflag) {
	    fprintf(stderr, "%s: -c and -f are incompatible\n", pname);
	    exit(1);
	}
	if (pflag) {
	    fprintf(stderr, "%s: -c and -p are incompatible\n", pname);
	    exit(1);
	}
	if (Pflag) {
	    fprintf(stderr, "%s: -c and -P are incompatible\n", pname);
	    exit(1);
	}
    }
    if (nflag & pflag & Pflag) {
	fprintf(stderr, "%s: -n, -p, and -P are incompatible\n", pname);
	exit(1);
    }
    if (Nflag) {
	if (nflag) {
	    fprintf(stderr, "%s: -n and -N are incompatible\n", pname);
	    exit(1);
	}
	else
	    nflag = Nflag;
    }
    if (smooth > 1) {
        if (cflag) {
	    fprintf(stderr, "%s: -c and -s or -S are incompatible\n", pname);
	    exit(1);
	}
        if (pflag) {
	    fprintf(stderr, "%s: -p and -s or -S are incompatible\n", pname);
	    exit(1);
	}
    }

    /* Make sure that len is a power of two. */
    if (len < 1) len = 1;
    if (len < LEN) {
	for (m = LEN; m >= len; m >>= 1)
	    ;
	m <<= 1;
    }
    else {
	for (m = LEN; m < len; m <<= 1)
	    ;
    }
    len = m;

    if ((c = (float *)calloc(len, sizeof(float))) == NULL) {
	fprintf(stderr, "%s: insufficient memory\n", pname);
	exit(2);
    }
	
    if (iflag) {		/* calculate and print inverse FFT */
	for (n = 0, rsum = 0.; n < len && fscanf(ifile, "%f", &c[n]) == 1; n++)
	    ;
	if (n == 0) {
	    fprintf(stderr, "%s: standard input is empty\n", pname);
	    exit(2);
	}
	ifft();
	exit(0);
    }

    else {			/* calculate and print forward FFT */
	if (nflag) {		/* process input in chunks */
	    float *s, *t;
	    int nf2 = nflag/2;

	    for (m = len; m >= nflag; m >>= 1)
		;
	    m <<= 1;		/* m is now the smallest power of 2 >= nflag */
	    if ((s = (float *)calloc(sizeof(float), m)) == NULL ||
		(t = (float *)calloc(sizeof(float), m/2)) == NULL) {
		fprintf(stderr, "%s: insufficient memory\n", pname);
		exit(2);
	    }
	    for (n = 0; n < nf2 && fscanf(ifile, "%f", &t[n]) == 1; n++)
		;
	    while (1) {
		if (zflag) {
		    for (n = 0, rsum = 0.; n < nf2; n++)
			rsum += c[n] = t[n];
		    for (i=0; n<nflag && fscanf(ifile,"%f",&t[i])==1; i++, n++)
			rsum += c[n] = t[i]; /* read input, accumulate sum */
		}
		else {
		    for (n = 0, rsum = 0.; n < nf2; n++)
			c[n] = t[n];
		    for (i=0; n<nflag && fscanf(ifile,"%f",&t[i])==1; i++, n++)
			c[n] = t[i];
		    ;
		}
		if (n < nflag) break;
		fft();
		if (Nflag)
		    fft_out();
		else if (Pflag)
		    for (i = 0; i < m; i++)
			s[i] += c[i]*c[i];
		else
		    for (i = 0; i < m; i++)
			s[i] += c[i];
		nchunks++;
	    }
	    if (nchunks < 1) {
		fprintf(stderr, "%s: input series is too short\n", pname);
		exit(2);
	    }
	    if (Nflag == 0 && Pflag)
		for (i = 0; i < m; i++)
		    c[i] = sqrt(s[i])/nchunks;
	    else
		for (i = 0; i < m; i++)
		    c[i] = s[i]/nchunks;
	}

	else {
	    read_input();
	    if (n == 0) {
		fprintf(stderr, "%s: standard input is empty\n", pname);
		exit(2);
	    }
	    fft();
	}

	if (!Nflag)
	    fft_out();

	exit(0);
    }
}

/* This function detrends (subtracts a least-squares fitted line from) a
   a sequence of n uniformily spaced ordinates supplied in c. */
detrend(c, n)
float *c;
int n;
{
    int i;
    double a, b = 0.0, tsqsum = 0.0, ysum = 0.0, t;

    for (i = 0; i < n; i++)
	ysum += c[i];
    for (i = 0; i < n; i++) {
	t = i - n/2 + 0.5;
	tsqsum += t*t;
	b += t*c[i];
    }
    b /= tsqsum;
    a = ysum/n - b*(n-1)/2.0;
    for (i = 0; i < n; i++)
	c[i] -= a + b*i;
    if (b < -0.04 || b > 0.04)
	fprintf(stderr,
		"%s: (warning) possibly significant trend in input series\n",
		pname);
}

read_input()
{
    if (zflag)
	for (n = 0, rsum = 0.; n < len && fscanf(ifile, "%f", &c[n]) == 1; n++)
	    rsum += c[n];	/* read input, accumulate sum */
    else
	for (n = 0, rsum = 0.; n < len && fscanf(ifile, "%f", &c[n]) == 1; n++)
	    ;
}

fft()		/* calculate forward FFT */
{
    int i;

    if (zflag) {		/* zero-mean the input array */
	rmean = rsum/n;
	for (i = 0; i < n; i++)
	    c[i] -= rmean;	
	if (zflag == 2)
	    detrend(c, n);
    }
    for (m = len; m >= n; m >>= 1)
	;
    m <<= 1;		/* m is now the smallest power of 2 >= n; this is the
			   length of the input series (including padding) */
    if (wflag)			/* apply the chosen windowing function */
	for (i = 0; i < m; i++)
	    c[i] *= (*window)(i, m);
    else
	wsum = m;
    norm = sqrt(2.0/(wsum*n));
    if (fflag) fstep = freq/(2.*m); /* note that fstep is actually half of
				       the frequency interval;  it is
				       multiplied by the doubled index i
				       to obtain the center frequency for
				       bin (i/2) */
    realft(c-1, m/2, 1);	/* perform the FFT;  see Numerical Recipes */
}

fft_out()	/* print the FFT */
{
    int i;

    c[m] = c[1];		/* unpack the output array */
    c[1] = c[m+1] = 0.;
    for (i = 0; i <= m; i += 2*decimation) {
	int j;
	double pow;

	if (fflag) printf("%g\t", i*fstep);
	if (cflag) printf("%g\t%g\n", c[i], c[i+1]);
	else {
	    for (j = 0, pow = 0.0; j < 2*smooth; j += 2)
	        pow += (c[i+j]*c[i+j] + c[i+j+1]*c[i+j+1])*norm*norm;
	    pow /= smooth/decimation;
	    if (Pflag) printf("%g", pow);
	    else printf("%g", sqrt(pow));
	    if (pflag) printf("\t%g", atan2(c[i+1], c[i]));
	    printf("\n");
	}
    }
}

ifft()		/* calculate and print inverse FFT */
{
    int i;

    n -= 2;
    c[1] = c[n];		/* repack IFFT input array */
    if (iflag < 0) {		/* convert polar form input to rectangular */
	for (i = 2; i < n; i += 2) {
	    float im;

	    im = c[i]*sin(c[i+1]);
	    c[i] *= cos(c[i+1]);
	    c[i+1] = im;
	}
    }
    realft(c-1, n/2, -1);
    if (iflag < 0) {
	norm = sqrt(2.0);
	for (i = 0; i < n; i++)
	    printf("%g\n", c[i]*norm);
    }
    else
	for (i = 0; i < n; i++)
	    printf("%g\n", c[i]/(n/2.0));
}

char *prog_name(s)
char *s;
{
    char *p = s + strlen(s);

#ifdef MSDOS
    while (p >= s && *p != '\\' && *p != ':') {
	if (*p == '.')
	    *p = '\0';		/* strip off extension */
	if ('A' <= *p && *p <= 'Z')
	    *p += 'a' - 'A';	/* convert to lower case */
	p--;
    }
#else
    while (p >= s && *p != '/')
	p--;
#endif
    return (p+1);
}

static char *help_strings[] = {
 "",
 NULL
};

void help()
{
    int i;

    (void)fprintf(stderr, help_strings[0], pname);
    for (i = 1; help_strings[i] != NULL; i++) {
	(void)fprintf(stderr, "%s\n", help_strings[i]);
	if (i % 23 == 0) {
	    char b[5];
	    (void)fprintf(stderr, "--More--");
	    (void)fgets(b, 5, stdin);
	    (void)fprintf(stderr, "\033[A\033[2K"); /* erase "--More--";
						       assumes ANSI terminal */
	}
    }
}

void realft(data,n,isign)
float data[];
int n,isign;
{
    int i, i1, i2, i3, i4, n2p3;
    float c1 = 0.5, c2, h1r, h1i, h2r, h2i;
    double wr, wi, wpr, wpi, wtemp, theta;
    void four1();

    theta = PI/(double) n;
    if (isign == 1) {
	c2 = -0.5;
	four1(data, n, 1);
    } 
    else {
	c2 = 0.5;
	theta = -theta;
    }
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0+wpr;
    wi = wpi;
    n2p3 = 2*n+3;
    for (i = 2; i <= n/2; i++) {
	i4 = 1 + (i3 = n2p3 - (i2 = 1 + ( i1 = i + i - 1)));
	h1r =  c1*(data[i1] + data[i3]);
	h1i =  c1*(data[i2] - data[i4]);
	h2r = -c2*(data[i2] + data[i4]);
	h2i =  c2*(data[i1] - data[i3]);
	data[i1] =  h1r + wr*h2r - wi*h2i;
	data[i2] =  h1i + wr*h2i + wi*h2r;
	data[i3] =  h1r - wr*h2r + wi*h2i;
	data[i4] = -h1i + wr*h2i + wi*h2r;
	wr = (wtemp = wr)*wpr - wi*wpi+wr;
	wi = wi*wpr + wtemp*wpi + wi;
    }
    if (isign == 1) {
	data[1] = (h1r = data[1]) + data[2];
	data[2] = h1r - data[2];
    } else {
	data[1] = c1*((h1r = data[1]) + data[2]);
	data[2] = c1*(h1r - data[2]);
	four1(data, n, -1);
    }
}

void four1(data, nn, isign)
float data[];
int nn, isign;
{
    int n, mmax, m, j, istep, i;
    double wtemp, wr, wpr, wpi, wi, theta;
    float tempr, tempi;
    
    n = nn << 1;
    j = 1;
    for (i = 1; i < n; i += 2) {
	if (j > i) {
	    tempr = data[j];     data[j] = data[i];     data[i] = tempr;
	    tempr = data[j+1]; data[j+1] = data[i+1]; data[i+1] = tempr;
	}
	m = n >> 1;
	while (m >= 2 && j > m) {
	    j -= m;
	    m >>= 1;
	}
	j += m;
    }
    mmax = 2;
    while (n > mmax) {
	istep = 2*mmax;
	theta = TWOPI/(isign*mmax);
	wtemp = sin(0.5*theta);
	wpr = -2.0*wtemp*wtemp;
	wpi = sin(theta);
	wr = 1.0;
	wi = 0.0;
	for (m = 1; m < mmax; m += 2) {
	    for (i = m; i <= n; i += istep) {
		j =i + mmax;
		tempr = wr*data[j]   - wi*data[j+1];
		tempi = wr*data[j+1] + wi*data[j];
		data[j]   = data[i]   - tempr;
		data[j+1] = data[i+1] - tempi;
		data[i] += tempr;
		data[i+1] += tempi;
	    }
	    wr = (wtemp = wr)*wpr - wi*wpi + wr;
	    wi = wi*wpr + wtemp*wpi + wi;
	}
	mmax = istep;
    }
}