/*
 * fft - DFT, IDFT using FFT
 *
 * 1970 Kouichi Hirabayashi
 * 1999-05-23 Add new windows
 * 2008-08-31 Add time/frequency field and step-like function correction.
 * 2009-01-10 Add chirp z transform.
 */

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

/* input option */
#define TIME		1

/* output option */
#define ABS		1
#define POWER		2
#define REAL		3
#define FILTER		4

/* window option */
#define HAMMING		1
#define HANNING		2
#define BLACKMAN	3
#define COS		4
#define RIESZ		5
#define FLAT_TOP	6

/* transfoem option */
#define DFT	0
#define IDFT	1
#define RAW	2

typedef struct dfts {
  int n;		/* # of data */
  int l;		/* 2^l == n */
  int no;		/* # of output data */
  int nx;		/* # of FFT data */
  int nxy;		/* data size */
  double *x, *y;	/* input data */
  double dt;		/* time domain sample interval */
  double df;		/* frequency domain sample interval */
  double u, v;		/* last data */
  double *wr, *wi;	/* n */
  double *vr, *vi;	/* nx */
  double *tr, *ti;	/* nx */
} dfts;

dfts dftd;

FILE *ifp, *ofp, *efopen();

double PI, PI2;
int cflg;	/* chirp z transform */
int sflg;	/* step-like function */
int topt;	/* transform option */
int iopt;	/* input option */
int oopt;	/* output option */
int wopt;	/* window option */

char buf[BUFSIZ], *cmd;

main(argc, argv) char **argv;
{
  char *s;

  cmd = argv[0];
  dftd.dt = dftd.df = 0.0;
  while (--argc > 0 && (*++argv)[0] == '-')
	for (s = argv[0]+1; *s != '\0'; s++)
		switch (*s) {
		case 'B':	/* Blackman-Harris window */
			wopt = BLACKMAN;
			break;
		case 'F':	/* flat-top window */
			wopt = FLAT_TOP;
			break;
		case 'M':	/* Hamming window */
			wopt = HAMMING;
			break;
		case 'N':	/* Hanning window */
			wopt = HANNING;
			break;
		case 'R':	/* Riesz window */
			wopt = RIESZ;
			break;
		case 'a':	/* abs, arg form output */
			oopt = ABS;
			break;
		case 'c':	/* chirp z transform */
			cflg = 1;
			break;
		case 'd':	/* DFT */
			topt = DFT;
			break;
		case 'f':	/* sampling interval (frequency) */
			if (*(s + 1) == '\0') {
				s = *++argv; argc--;
			}
			else
				s++;
			sscanf(s, "%le", &dftd.df);
			s += strlen(s) - 1;
			break;
		case 'i':	/* IDFT */
			topt = IDFT;
			break;
		case 'n':	/* filter only */
			oopt = FILTER;
			break;
		case 'p':	/* power spectram (dB) output */
			oopt = POWER;
			break;
		case 'r':	/* FFT only */
			topt = RAW;
			break;
		case 's':	/* Step-Like function */
			sflg++;
			break;
		case 't':	/* sampling interval (time) */
			if (*(s + 1) == '\0') {
				s = *++argv; argc--;
			}
			else
				s++;
			sscanf(s, "%le", &dftd.dt);
			s += strlen(s) - 1;
			break;
		case 'x':
			oopt = REAL;
			break;
		default:
			Usage();
		}
  ifp = (argc > 0 || (argc > 0 && strcmp(argv[0],  "-"))) ?
	efopen(argv[0], "r") : stdin;
  ofp = (argc > 1) ? efopen(argv[1], "w") : stdout;
  if (topt == RAW)
	cflg = 0;
  PI = atan2(0.0, -1.0);
  PI2 = PI * 2.0;

  dft_init(&dftd, ifp);
  if (wopt)
	filter(&dftd);

  if (cflg) {
	 ChirpZ(&dftd);
  }
  else
	Dft(&dftd);

  output(&dftd);
  dft_clean(&dftd);
  return 0;
}

/* DFT using FFT */

Dft(dfts *p)
{
  int i;
  double u, v;

  fft(p->x, p->y, p->l);
  if (sflg)
	step(p);
  if (topt == DFT) {
	u = (p->dt > 0.0) ? p->dt : 1.0 / (double)p->nx;
	v = -u;
	for (i = 0; i < p->nx; i++) {
		p->x[i] *= u;
		p->y[i] *= v;
	}
  }
  else if (topt == IDFT) {	/* IDFT */
	u = (p->df > 0) ? p->df : 1.0;
	v = -u;
	if (p->df > 0) {
		for (i = 0; i < p->nx; i++) {
			p->x[i] *= u;
			p->y[i] *= v;
		}
	}
  }
  return 0;
}

output(dfts *p)
{
  int i;
  double u, v, w, z;

  if (topt == RAW) {
	if (oopt == REAL)
		for (i = 0; i < p->no; i++)
			fprintf(ofp, "%g\n", p->x[i]);

	else
		for (i = 0; i < p->no; i++)
			fprintf(ofp, "%g %g\n", p->x[i], p->y[i]);
  }
  else if (topt == DFT) {
	w = (p->dt > 0) ? 1.0 / (double)p->n / p->dt : 1.0 / (double)p->n;
	if (oopt == POWER) {
		for (i = 0; i < p->no; i++) {
			u = p->x[i]; v = p->y[i];
			z = u * u + v * v;
			z = 20.0 * log10(z);
			if (p->dt > 0.0)
				fprintf(ofp, "%g ", w * (double)i);
			fprintf(ofp, "%g\n", z);
		}
	}
	else if (oopt == ABS) {
		for (i = 0; i < p->no; i++) {
			u = p->x[i]; v = p->y[i];
			if (u != 0.0 || v != 0.0) {
				if (p->dt > 0.0)
					fprintf(ofp, "%g ", w * (double)i);
				fprintf(ofp, "%g %g",
					sqrt(u * u + v * v), atan2(u, v));
			}
			else {
				if (p->dt > 0.0)
					fprintf(ofp, "%g ", w * (double)i);
				fprintf(ofp, "%g %g",
					sqrt(u * u + v * v), 0.0);
			}
		}
	}
	else {
		for (i = 0; i < p->no; i++) {
			if (p->dt > 0.0)
				fprintf(ofp, "%g ", w * (double)i);
			fprintf(ofp, "%g %g\n", p->x[i], p->y[i]);
		}
	}
  }
  else if (topt == IDFT) {
	w = (p->df > 0.0) ? 1.0 / (double)p->n / p->df : 1.0 / (double)p->n;
	for (i = 0; i < p->no; i++) {
		if (p->df > 0)
			fprintf(ofp, "%g ", w * (double)i);
		if (oopt == REAL)
			fprintf(ofp, "%g\n", p->x[i]);
		else
			fprintf(ofp, "%g %g\n", p->x[i], p->y[i]);
	}
  }
  return 0;
}

dft_init(dfts *p, FILE *fp)
{
  int i, j, k, j0, n;
  double start, prev, dmin, dmax;
  char *s;
  void *emalloc(), *erealloc();
  double u, v, w, z, u0, v0;

  n = 0;
  while ((s = fgets(buf, sizeof(buf), ifp)) != NULL) {
	if (*buf == '#')
		continue;
	else if (strncmp(buf, ".T", 2) == 0) {
		sscanf(buf, "%*s %le %d", &p->dt, &p->n);
		topt = DFT;
	}
	else if (strncmp(buf, ".F", 2) == 0) {
		sscanf(buf, "%*s %le %d", &p->df, &p->n);
		topt = IDFT;
	}
	else
		break;
  }
  if (s == NULL)
	exit(1);

  p->nxy = 1024;
  p->x = (double *)erealloc((double *) p->x,
		(unsigned)(p->nxy * sizeof(double)));
  p->y = (double *)erealloc((double *) p->y,
		(unsigned)(p->nxy * sizeof(double)));
  dmin = dmax = 0.0;
  for (i = 0; i < p->nxy; ) {
	if (*buf == '#')
		continue;
	u = v = w = 0.0;
	j = sscanf(buf, "%le %le %le", &u, &v, &w);
	if (i == 0) {	/* the first data line */
		u0 = u;
		v0 = v;
		if (j == 2 && topt == DFT)
			iopt = TIME;
	}
	if (j == 3 || iopt == TIME) {	/* time series */
		if (i == 0)
			start = prev = u;
		else if (i == 1) 
			dmin = dmax = u - start, prev = u;
		else {
			z = u - prev;
			prev = u;
			if (z < dmin)
				dmin = z;
			if (dmax < z)
				dmax = z;
		}
		p->x[i] = v;
		p->y[i] = w;
	}
	else {
		p->x[i] = u;
		p->y[i] = v;
	}
	if (i > 0 && j != j0)
		fprintf(stderr, "Warning: # of field error: %s", buf);
	j0 = j;
	i++;
	if (fgets(buf, BUFSIZ, ifp) == NULL)
		break;
	if (i == p->nxy - 1) {	/* memory allocation */
		p->nxy <<= 1;
		p->x = (double *)erealloc((double *) p->x,
			(unsigned)(p->nxy * sizeof(double)));
		p->y = (double *)erealloc((double *) p->y,
			(unsigned)(p->nxy * sizeof(double)));
	}
  }
  if (dmin + dmax > 0) {
	if (fabs(dmax - dmin) / (dmin + dmax) > 0.01)
		fprintf(stderr,
			"Warning: abnormal sampling interval: %g, %g\n",
				dmin, dmax);
  }
  if (j < 3 && iopt != TIME)
	w = 1.0;
  else
	w = (u - u0) / (double)(i - 1);
  if (p->dt == 0.0 && topt == DFT)
	p->dt = w;
  else if (p->df == 0.0 && topt == IDFT)
	p->df = w;

  p->n = i;	/* # of data */
  if (!cflg) {	/* FFT */
	p->l = (int)(log((double)i) / log(2.0) + 0.001);
	k = pow(2.0, (double)p->l);
	if ((topt == DFT && k < i) || (topt == IDFT && k < i - 1))
		p->l++;
	p->nx = (int)(pow(2.0, (double)p->l));
	/* fill data */
	if (sflg) {	/* last data */
		p->u = p->x[i - 1];
		p->v = p->y[i - 1];
	}
	else
		p->u = p->v = 0.0;
	/* stretch or zero fill */
	for ( ; i < p->nx; i++) {
		p->x[i] = p->u;
		p->y[i] = p->v;
	}
	p->no = (topt == RAW) ? p->nx : p->nx / 2 + 1;
  }
  else if (topt == DFT) {	/* chirp z transform */
	p->nx = p->n;
	p->no = p->n / 2;
  }

  /* add latter half (minus frequency) */
  if (topt == IDFT) {
	k = p->nx;
	if (p->nxy < 2 * k) {
		p->nxy = 2 * k;
		p->x = (double *)erealloc((double *)p->x,
			(unsigned)(p->nxy * sizeof(double)));
		p->y = (double *)erealloc((double *)p->y,
			(unsigned)(p->nxy * sizeof(double)));
	}
	if (cflg) {
		k = p->n - 1;
		n = k << 1;
		p->n = n;
		p->no = n;
		/* add latter half (minus frequency) */
		for (i = 1; i < k; i++) {
			p->x[n - i] = p->x[i];
			p->y[n - i] = -p->y[i];
		}
		p->y[k] = 0;
		return 0;
	}
	n = 2 * p->nx;
	/* add latter half (minus frequency) */
	for (i = 1; i < k; i++) {
		p->x[n - i] = p->x[i];
		p->y[n - i] = -p->y[i];
	}
	p->y[k] = 0;
	p->nx = n;
	p->l++;
	p->no = n;
	p->n = n;
  }
  return 0;
}

dft_clean(dfts *p)
{
  xfree(&p->x);
  xfree(&p->y);
  xfree(&p->wr);
  xfree(&p->vr);
  xfree(&p->tr);
  p->n = p->no = p->nx = p->l = 0;
  return 0;
}

xfree(p) void **p;
{
  if (*p != NULL) {
	free(*p);
	*p = NULL;
  }
}

filter(dfts *p)
{
  int i;
  double w;

  if (wopt == BLACKMAN)
	blackman(p->x, p->y, p->n);
  else if (wopt == HANNING)
	hanning(p->x, p->y, p->n);
  else if (wopt == HAMMING)
	hamming(p->x, p->y, p->n);
  else if (wopt == COS)
	half_cos(p->x, p->y, p->n);
  else if (wopt == RIESZ)
	riesz(p->x, p->y, p->n);
  else if (wopt == FLAT_TOP)
	flat_top(p->x, p->y, p->n);
  else
	fprintf(stderr, "Undefined windows (%d)\n", wopt);
  if (oopt == FILTER) {	/* filter only */
	if (topt == DFT)
		w = (p->dt > 0) ? 1.0 / (double)p->n / p->dt : 1.0;
	else
		w = (p->df > 0) ? 1.0 / (double)p->n / p->df : 1.0;
	for (i = 0; i < p->n; i++) {
		if (iopt == TIME)
			fprintf(ofp, "%g %g\n", w * (double)i, p->x[i]);
		else
			fprintf(ofp, "%g %g %g\n",
				w * (double)i, p->x[i], p->y[i]);
	}
	dft_clean(&dftd);
	exit(0);
  }
  return 0;
}

/*
 * Step-Like functions correction
 *
 * Add (u - u0) / (1 - exp(-j*2*π/N) to DFT result.
 */
step(dfts *p)
{
  int i;
  double a, u, v, w;

  u = cflg ? -p->u : p->u;
  v = PI2 / (double)p->n;
  for (i = 1; i < p->n; i++) {
	w = (double)i * v;
	p->x[i] += 0.5 * p->u;
	p->y[i] += 0.5 * u * sin(w) / (1.0 - cos(w));
  }
}

void *emalloc(i) unsigned i;
{
  char *p;

  if ((p = malloc(i)) == NULL)
	error("memory over");
  return p;
}

void *erealloc(p, i) char *p; unsigned i;
{
  char *q;

  if ((q = realloc(p, i)) == NULL)
	error("memory over", (char *) 0);
  return q;
}

FILE *efopen(file, mode) char *file, *mode;
{
  FILE *fp;

  if ((fp = fopen(file, mode)) == NULL)
	error("can't open %s", file);
  return fp;
}

/*
 * fast fourier transform
 * x: real part
 * y: imaginal part
 * l: # of data = 2 ^ l
 */
fft(x, y, l) double *x, *y; int l;
{
  int i, j, k, n, lm, lmx, l0, li, lix, nm1, nv2;
  double c, s, t1, t2, arg, scl, sin(), cos();

  n = 1 << l;
  lmx = n;
  scl = PI2 / n;
  for (l0 = 0; l0 < l; scl *= 2.0, l0++) {
	lix = lmx, lmx /= 2, arg = 0.0;
	for (lm = 0; lm < lmx; lm++) {
		c = cos(arg); s = -sin(arg); arg += scl;
		for (li = lix; li <= n; li += lix) {
			i = li - lix + lm;
			j = i + lmx;
			t1 = *(x + i) - *(x + j);
			t2 = *(y + i) - *(y + j);
			*(x + i) += *(x + j);
			*(y + i) += *(y + j);
			*(x + j) = c * t1 + s * t2;
			*(y + j) = c * t2 - s * t1;
		}
	}
  }
  nv2 = n / 2; nm1 = n - 1; j = 0;
  for (i = 0; i < nm1; j += k, i++) {
	if (i < j) {
		t1 = *(x + j); t2 = *(y + j);
		*(x + j) = *(x + i); *(y + j) = *(y + i);
		*(x + i) = t1; *(y + i) = t2;
	}
	for (k = nv2; k <= j; j -= k, k /= 2)
		;
  }
}

/*
 * Windows
 */

/* Hanning window (power reduction 4.26 dB) */

hanning(x, y, n) double *x, *y; int n;
{
  int i;
  double u, v;

  u = PI2 / (double) n;
  for (i = 0; i < n; i++) {
	v = 0.50 * (1.0 + cos((double)i * u));
	*(x + i) *= v;
	*(y + i) *= v;
  }
}

/* Hamming window (power reduction 4.0 dB) */

hamming(x, y, n) double *x, *y; int n;
{
  int i;
  double u, v;

  u = PI2 / (double) n;
  for (i = 0; i < n; i++) {
	v = 0.54 + 0.46 * cos((double)i * u);
	*(x + i) *= v;
	*(y + i) *= v;
  }
}

/* Blackman-Harris window (power reduction 5.14 dB) */

blackman(x, y, n) double *x, *y; int n;
{
  int i;
  double u1, u2, v;

  u1 = PI2 / (double) n;
  u2 = 2.0 * u1;
  for (i = 0; i < n; i++) {
	v = 0.423 + 0.498 * cos((double)i * u1)
		+ 0.0792 * cos((double)i * u2);
	*(x + i) *= v;
	*(y + i) *= v;
  }
}

/* Halh cosine window (power reduction 3.3 dB) */

half_cos(x, y, n) double *x, *y; int n;
{
  int i;
  double u, v;

  u = PI / (double) n;
  for (i = 0; i < n; i++) {
	v = sin((double)i * u);
	*(x + i) *= v;
	*(y + i) *= v;
  }
}

/* Riesz window (power reduction 2.73 dB) */

riesz(x, y, n) double *x, *y; int n;
{
  int i;
  double u, v;

  u = 1.0 / (double)n / (double)n;
  for (i = 0; i < n; i++) {
	v = 4.0 * (double)i * (double)(n - i) * u;
	*(x + i) *= v;
	*(y + i) *= v;
  }
}

/* flat-top windows (power reduction 7.0 dB) */

flat_top(x, y, n) double *x, *y; int n;
{
  int i;
  double u1, u2, v, v1;

  u1 = PI2/(double)n;
  u2 = 2.0 * u1;
  for (i = 0; i < n; i++) {
	v = 0.54 - 0.46 * cos((double)i * u1);
	if (i > 0)
		v *= sin((double)i * u2) / (double)i / u2;
	*(x + i) *= v;
	*(y + i) *= v;
  }
}

/* chirp z transform */

ChirpZ(dfts *p)
{
  int i, n2;
  double u, v, w;

  if (czt_init(p) != 0)
	return -1;
  if (topt == DFT) {
	p->no++;
	u = p->x[p->n - 1];
	p->u = p->x[p->n - 1];
	p->v = p->y[p->n - 1];
	czt(&dftd, 0);
	if (sflg)
		step(p);
	u = (p->dt > 0.0) ? p->dt : 1.0 / (double)p->n;
	v = u;
	w = (p->dt > 0.0) ? 1.0 / (double)p->n / p->dt : 1.0 / (double)p->n;
	for (i = 0; i <= p->n; i++) {
		p->x[i] *= u;
		p->y[i] *= v;
	}
#if 0
	n2 = p->n / 2;
	fprintf(stderr, "n2(%d)\n", n2);
	for (i = 0; i <= n2; i++) {
		if (p->dt > 0.0)
			fprintf(ofp, "%g ", w * (double)i);
		fprintf(ofp, "%g %g\n", p->x[i], p->y[i]);
	}
	dft_clean(p);
	exit(0);
#endif
  }
  else {
	czt(&dftd, 1);
	if (p->df > 0.0)
		u = p->df * (double)p->n;
	else
		u = 1.0;
	for (i = 0; i < p->n; i++) {
		p->x[i] *= u;
		p->y[i] *= -u;
	}
#if 0
	w = (p->df > 0.0) ? 1.0 / u : 1.0 / (double)p->n;
	for (i = 0; i < p->no; i++) {
		if (p->df > 0)
			fprintf(ofp, "%g ", w * (double)i);
		fprintf(ofp, "%g %g\n", p->x[i], p->y[i]);
	}
	dft_clean(p);
	exit(0);
#endif
  }
  return 0;
}

/*
 * 重みデータ・インパルス応答データを作る
 */
//static void make_cztdata(int n, int no, int nx, double wr[], double wi[],
static void make_cztdata(int n, int nx, double wr[], double wi[],
                         double vr[], double vi[])
{
  double r, d = PI / n;
  int i, j, n2 = n * 2;

  /*
   * W^nk
   * W = exp(-j*2π/n)
   * nk = (n^2 + k^2 - (k - n)^2) / 2
   *
   * w(k) = W^(k^2/2)
   * v(k) = W^(k^2/2)
   */
  for (i = j = 0; i < n; i++, j = (j + (2 * i - 1)) % n2) {
	r = d * j;
	wr[i] =  cos(r);			/*  cos(i*i*PI/n) */
	wi[i] = -sin(r);			/* -sin(i*i*PI/n) */
  }

  vr[0] = 1.0;
  vi[0] = 0.0;
  for (i = 1; i < n; i++) {
	vr[i] = vr[nx - i] =  wr[i];	/* nx * cos(i*i*PI/n) */
	vi[i] = vi[nx - i] = -wi[i];	/* nx * sin(i*i*PI/n) */
  }
//  for (i = no, j = nx - n; i <= j; i++, j--) {
  for (i = n, j = nx - n; i <= j; i++, j--) {
	vr[j] = vr[i] = 0;			/* 残りの拡張部には 0 を */
	vi[j] = vi[i] = 0;
  }
}

/*
 * CZT計算用構造体に対し、標本数 n, 出力データ数 no 用の数表データを
 * 作成する。
 *
 * cztp	= CZT計算用構造体へのポインタ
 * n		= 標本点の数
 * no		= 出力するデータの数
 * return	= 0:正常終了 1:nが無効な数 2:メモリ不足
 */
int czt_init(dfts *p)
{
  int i;

  if (p->n <= 1)
	return -1;
  if (p->no <= 1 || p->n < p->no)
	p->no = p->n;

  p->nx = p->n + p->no;
  /* p->nx を2の整数乗まで拡張する */
  for (i = 1; i < p->nx; i *= 2)
	;
  p->nx = i;

  p->wr = malloc(2 * p->n  * sizeof(double));
  p->vr = malloc(2 * p->nx * sizeof(double));
  p->tr = malloc(2 * p->nx * sizeof(double));
  if (p->wr == NULL || p->vr == NULL || p->tr == NULL) {
	dft_clean(p);
	return 2;
  }
  p->wi = p->wr + p->n;
  p->vi = p->vr + p->nx;
  p->ti = p->tr + p->nx;

  /* V(n) */
//  make_cztdata(n, no, nx, p->wr, p->wi, p->vr, p->vi);
  make_cztdata(p->n, p->nx, p->wr, p->wi, p->vr, p->vi);
  /* FFT */
  fft(p->vr, p->vi, (int)(log(p->nx) / log(2.0) + 0.001));
  return 0;
}

/*
 * CZT (Chirp z-Tramsform) アルゴリズムによる高速 DFT.
 * re[] が実部, im[] が虚部. 結果は re[],im[] に上書きされる.
 * inv!=0 (=TRUE) なら逆変換を行う. cztp には, 計算用データが
 * 入っている構造体を指定する.
 */
//void czt(dfts *p, int inv, double re[], double im[])
//czt(dfts *p, int inv, double re[], double im[])
czt(dfts *p, int inv)
{
  double xr, xi, yr, yi, t, *re, *im;
  int i, n, no, nx;

  re = p->x;
  im = p->y;
  n = p->n;
  no = p->no;
  nx = p->nx;

  /* Y(k) = g(k) * W^(k^2/2) */
  for (i = 0; i < n; i++) {		/* 入力の重みデータ乗算 */
	yr = p->wr[i];
	yi = p->wi[i];
	if (inv) yi = -yi;
	p->tr[i] = re[i] * yr - im[i] * yi;
	p->ti[i] = -(im[i] * yr + re[i] * yi);
  }
  for (; i < nx; i++) {			/* 残りの拡張部には 0 を */
	p->tr[i] = 0;
	p->ti[i] = 0;
  }
  /* FFT */
  fft(p->tr, p->ti, (int)(log(p->nx) / log(2.0) + 0.001));

  /* Y(k)*V(k) */
  for (i = 0; i < p->nx; i++) {		/* コンボリューション演算 */
	xr = p->tr[i];
	xi = -p->ti[i];
	yr = p->vr[i];
	yi = p->vi[i];
	if (inv) yi = -yi;
	p->tr[i] = xr * yr - xi * yi;
	p->ti[i] = xi * yr + xr * yi;
  }
  /* IFFT */
  fft(p->tr, p->ti, (int)(log(p->nx) / log(2.0) + 0.001));
  for (i = 0; i < nx; i++) {
	p->tr[i] = p->tr[i] / nx;
	p->ti[i] = p->ti[i] / nx;
  }

  /* IFFT(Y(k)*V(k)) * W^(n^2/2) */
  for (i = 0; i < p->no; i++) {		/* 出力の重みデータ乗算 */
	yr = p->wr[i];
	yi = p->wi[i];
	if (inv) yi = -yi;
	re[i] = p->tr[i] * yr - p->ti[i] * yi;
	im[i] = p->ti[i] * yr + p->tr[i] * yi;
  }

  if (inv) {				/* 逆変換ならnで割る */
	t = 1.0 / n;			/* 逆数をかける(除算は遅いので) */
	for (i = 0; i < p->no; i++) {
		re[i] *= t;
		im[i] *= t;
	}
  }
  return 0;
}

error(s, t)
char *s, *t;
{
  fprintf(stderr, "%s: ", cmd);
  fprintf(stderr, s, t);
  fprintf(stderr, "\n");
  exit(1);
}

Usage()
{
  static char *t[] = {
"コマンド	fft - FFT による離散フーリエ変換 (DFT)",
"構文		fft [オプション] [入力ファイル [出力ファイル]]",
"オプション	-B	入力データに Blackman-Harris フィルタを掛けてから変換",
"		-F	入力データに flat-top フィルタを掛けてから変換",
"		-M	入力データに Hammin フィルタを掛けてから変換",
"		-N	入力データに Hanning フィルタを掛けてから変換",
"		-R	入力データに Riesz フィルタを掛けてから変換",
"		-a	出力データを絶対値と偏角で表示",
"		-c	(FFT でなく) chirp z 変換を使います",
"		-f F	周波数軸の標本化間隔 F (Hz) を指定",
"		-i	フーリエ逆変換(IDFT)を実行",
"		-n	フィルタリングのみ実行",
"		-p	出力データをパワー(dB)で表示",
"		-r	FFT 演算のみ実行",
"		-s	ステップ状関数補正を実行",
"		-t T	時間軸の標本化間隔 T (s) を指定",
"		-x	-i を指定したとき、実部だけを表示",
"機能	-r, -i オプションにより次の３つの演算のいずれかを選択します。",
"  -r:     FFT 演算を実行",
"  -i:     離散フーリエ逆変換を実行",
"  その他: 離散フーリエ変換を実行",
"入力データの形式は次のとおりです。",
"  FFT 演算のとき: x y",
"  離散フーリエ変換のとき: t x [y]",
"  離散フーリエ変換のとき: f x y",
"  ここに、",
"	t = 時刻",
"	f = 周波数",
"	x = 関数の実部",
"	y = 関数の虚部 (省略すると 0)",
"	j = sqrt(-1)",
"ただし、-t, -f オプションか .T, .F データでサンプリング間隔を指定すれば、デー",
"タの t, f を省略して、x, y だけにすることもできます。この他下記形式のコメント",
"行とオプション行が使えます。",
"  [# コメント]",
"  [.T 時間軸の標本化間隔(s)]",
"  [.F 周波数軸の標本化間隔(Hz)]",
"時刻や周波数のサンプリング間隔は同じでなければなりません。また、周波数データ",
"の最初は", "0 Hz でなければなりません。データ数 N が２の巾乗でない場合は、-s",
"オプションが指定されていてば最後のデータを延長し、そうでなければ 0 を補って",
"２の巾乗にします。",
"出力データも同じ形式ですが、出力オプション -a, -p により、極座標形式にしたり、",
"パワースペクトル表示にすることができます。",
"演算の内容は、まず -r を指定した場合は FFT 演算",
"               N-1",
"  H(n) = (1/N)*Σ(x(k) + sqrt(-1)*y(k)]*exp(-sqrt(-1)*2*π*n*k/N)",
"               k=0",
"                             (n = 0, 1, 2, .. N-1)",
"を実行します。いずれの場合も、入力フィルタは実行されます。FFT 演算では、",
" 1) FFT 出力の虚部は連続フーリエ変換の虚部と符号が逆",
" 2) FFT 出力の前半 N/2+1 個が連続フーリエ変換のプラス側、後半がマイナス側",
"であることに注意してください。",
"以下、-r を指定しない場合になりますが、-i を指定しないときは、入力データに対",
"して",
"             N-1",
"  H(k*F) = T*Σ(x(k) + sqrt(-1)*y(k)]*exp(-sqrt(-1)*2*π*n*k/N)",
"             k=0",
"        (T = 時間軸サプリング間隔, F = 1/(N*T), n = 0, 1, 2, .. N-1)",
"",
"の離散フーリエ変換(DFT)演算を行い、その結果を指定された出力形式で表示します。",
"T を指定しないときは、T = 1/n になります。出力データはその虚部の符号を反転し、",
"正の周波数領域 N/2+1 個のデータだけを出力します。",
"-i を指定すると、入力データに対して、離散フーリエ逆変換(IDFT)",
"             N-1",
"  h(k*T) = F*Σ(x(k) + sqrt(-1)*y(k)*exp(sqrt(-1)*2*π*n*k/N)",
"             k=0",
"        (F = 周波数軸サンプリング間隔, T = 1/(N*F), n = 0, 1, 2, .. N-1)",
"を実行します。F を指定しないときは T の指定があれば F = 1/(N*T)、なければ、",
"F = 1、 になります。入力データは周波数軸のプラス側だけが指定されたものと見倣",
"し、実部が偶関数、虚部が奇関数になるようにデータを補います。",
"Hanning, Hamming, Blackman-Harris, Riesz, flat-top 等の入力フィルタは変換前に",
"実行されます。",
  NULL
  };
  int i;

  for (i = 0; t[i] != NULL; i++)
	fprintf(stdout, "%s\n", t[i]);
  exit(1);
}