/*--------------------------------------------------------------------
 *	$Id: spectrum1d.c,v 1.3.4.3 2002/02/27 17:58:55 pwessel Exp $
 *
 *	Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith
 *	See COPYING file for copying and redistribution conditions.
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 *
 *	Compute auto and cross spectra using Welch's method of
 *	multiple overlapped windows.  Find 1 standard error bars
 *	following expressions in Bendat & Piersol.
 *
 *	-D<dt>	set delta_time, the sampling of the timeseries.
 *		[Default <dt> = 1.0]
 *
 *	-N<name_stem>	name stem for filenames.  Files will be
 *		created called <name_stem>.xpower, etc.
 *		[Default <name_stem> = "spectrum"]
 *
 *	-S<segment_size>	give an integer radix 2 window width.
 *		Total timeseries will be split into pieces of
 *		length <segment_size>.  Std errors in spectra are
 *		approximately 1/sqrt(n_data/segment_size).
 *
 *	-V	Verbose operation; write info to stderr.
 *		[Default is silent]
 *
 *	-W	write Wavelength in col 1 instead of frequency.
 *		[Default writes frequency in cycles/dt_units]
 *
 *	-C<output flags>	input has 2 cols, X(t) Y(t); do cross-spectra.
 *		[Default is one column; do X power spectrum only]
 *		Optional string of 0 to 8 output flags { x y c n p a g o }
 *		0 or 8 produces all cross-spectra in output files.  x = x-power
 *		spectrum, y = y-power spectrum, c = coherent power spectrum,
 *		n = noise power spectrum, p = phase spectrum, a = admittance
 *		power spectrum, g = gain spectrum, o = squared coherency.
 *
 *	-b	input and/or output data are in binary form
 *
 *	Author:		W. H. F. Smith
 *	Date:		11 April 1991-2000
 *	Revised:	5 June 1991-2000 to add W and N options in prep for
 *			GMT v2.0 release.
 *			PW: Upgrade to GMT 3.1 w/ -b
 *			BCH-J: Upgraded -C: optional output flags
 *	References:	Julius S. Bendat & Allan G. Piersol,
 *			"Random Data", 2nd revised edition, 566pp.,
 *			1986, John Wiley & Sons, New York.
 *
 *			Peter D. Welch, "The use of Fast Fourier
 *			Transform for the estimation of power spectra:
 *			a method based on time averaging over short,
 *			modified periodograms", IEEE Transactions on
 *			Audio and Electroacoustics, Vol AU-15, No 2,
 *			June, 1967.
 *  Version:	3.4.1
 */

#include "gmt.h"

#define N_OUTPUT_CHOICES 8

float	*datac, *x, *y;

struct	SPEC {
		double	xpow;	/* PSD in X(t)  */
		double	ypow;	/* PSD in Y(t)  */
		double	gain;	/* Amplitude increase X->Y in Optimal Response Function  */
		double	phase;	/* Phase of Y w.r.t. X in Optimal Response Function  */
		double	coh;	/* (squared) Coherence of Y and X; SNR of Y = coh/(1-coh)  */
		double	radmit;	/* Real part of Admittance; used e.g., for gravity/topography  */
}	*spec;

double	dt = 1.0, x_variance, y_variance, d_n_windows, y_pow;

int	y_given = FALSE, write_wavelength = FALSE;
int	n_spec,  window = 0, window_2, n_data;

char	namestem[80], format[BUFSIZ];

FILE	*fp = NULL;

void	alloc_arrays(void);
void	read_data(void);
void	compute_spectra(void);
void	detrend_and_hanning(void);
void	write_output(char *output_choice, int n_outputs);
void	free_space(void);

main (int argc, char **argv)
{

	int	i, j, n_outputs, window_test, error = FALSE;
	char	output_choice[N_OUTPUT_CHOICES];

	argc = GMT_begin (argc, argv);

	n_outputs = 0;	
	for (i=0; i < N_OUTPUT_CHOICES; i++) output_choice[i] = 0;
	sprintf (namestem, "spectrum");

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
		
				/* Common parameters */
			
				case 'H':
				case 'V':
				case '\0':
					error += GMT_get_common_args (argv[i], 0, 0, 0, 0);
					break;
				
				/* Supplemental parameters */
				
				case 'b':
					error += GMT_io_selection (&argv[i][2]);
					break;
				case 'C':
					y_given = TRUE;
					j = 2;
				       while(argv[i][j]) {
					       switch (argv[i][j]) {
						       case 'x':	/* x power spectrum */
							       output_choice[j-2] = 'x';
							       break;
                                                       case 'y':	/* y power spectrum */
							       output_choice[j-2] = 'y';
                                                               break;
                                                       case 'c':	/* coherenct power spectrum */
							       output_choice[j-2] = 'c';
                                                               break;
                                                       case 'n':	/* noise power spectrum */
							       output_choice[j-2] = 'n';
                                                               break;
                                                       case 'p':	/* phase spectrum */
							       output_choice[j-2] = 'p';
                                                               break;
                                                       case 'a':	/* admittance spectrum */
							       output_choice[j-2] = 'a';
                                                               break;
                                                       case 'g':	/* gain spectrum */
							       output_choice[j-2] = 'g';
                                                               break;
                                                       case 'o':	/* (squared) coherency spectrum */
							       output_choice[j-2] = 'o';
                                                               break;
						       default:
							       error = TRUE;
							       fprintf (stderr, "%s: GMT SYNTAX ERROR -C option.  Unrecognized output choice %c\n", GMT_program, argv[i][j]);
					      }
					      j++;
					      n_outputs++;
					}
					if (j==2) { /* just -C: output all */
					      output_choice[0]='x';
					      output_choice[1]='y';
					      output_choice[2]='c';
					      output_choice[3]='n';
					      output_choice[4]='p';
					      output_choice[5]='a';
					      output_choice[6]='g';
					      output_choice[7]='o';
					      n_outputs=8;
				        }
					break;
				case 'D':
					dt = atof(&argv[i][2]);
					break;
				case 'N':
					strcpy(namestem, &argv[i][2]);
					if (!namestem[0]) error = TRUE;
					break;
				case 'S':
					window = atoi(&argv[i][2]);
					window_test = 2;
					while (window_test < window) {
						window_test += window_test;
					}
					break;
				case 'W':
					write_wavelength = TRUE;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else {
			if ((fp = GMT_fopen(argv[i], GMT_io.r_mode)) == NULL) {
				fprintf(stderr,"%s:  Cannot open r %s\n", GMT_program, argv[i]);
				error = TRUE;
			}
		}
	}

	if (argc == 1 || GMT_quick) {	/* Display usage */
		fprintf (stderr, "spectrum1d %s - compute auto- [and cross- ] spectra from one [or two] timeseries\n\n", GMT_VERSION);
		fprintf(stderr,"usage:  spectrum1d -S<segment_size> [-C[<xycnpago>]] [-D<dt>] [-H[<nrec>]] [-N<name_stem>] [-V] [-W] [-bi[s][<n>]] [-bo[s][<n>]]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf(stderr,"\t-S Use data subsets of <segment_size> elements.\n");
		fprintf(stderr,"\t   <segment_size> must be radix 2;\n");
		fprintf(stderr,"\t   std. err. = 1/sqrt(n_data/segment_size).\n");
		fprintf(stderr,"\tOptions:\n");
		fprintf(stderr,"\t-C[<xycnpago>] 2 column X(t),Y(t) input; estimate Cross-spectra\n\t   [Default 1 col, X power only].\n");
		fprintf(stderr,"\t   Optionally specify cross-spectra output(s)  [Default is all].\n");
		fprintf(stderr,"\t   x = xpower, y = ypower, c = coherent power, n = noise power\n");
		fprintf(stderr,"\t   p = phase, a = admittance, g = gain, o = squared coherency.\n\n");
		fprintf(stderr,"\t-D set delta_time sampling interval of data [Default = 1.0].\n");
		GMT_explain_option ('H');
		fprintf(stderr,"\t-N supply name stem for files [Default = 'spectrum'].\n");
		fprintf(stderr,"\t   Output files will be named <name_stem>.xpower, etc.\n");
		GMT_explain_option ('V');
		fprintf(stderr,"\t-W write Wavelength of spectral estimate in col 1 [Default = frequency].\n");
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf(stderr,"\t   Default is 2 input columns.\n");
		GMT_explain_option ('o');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}

	if (window <= 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: segment size must be positive\n", GMT_program);
		error++;
	}
	if (window_test != window) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: Segment size not radix 2.  Try %d or %d\n", GMT_program,
			(window_test/2), window_test);
		error++;
	}
	if (dt <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -D option: Sampling interval must be positive\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[0] && gmtdefs.io_header) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
        if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 2;
        if (GMT_io.binary[0] && GMT_io.ncol[0] < (y_given + 1)) {
                fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data must have at least %d columns\n", GMT_program, y_given + 1);
		error++;
	}

	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands */
	
	if (GMT_io.binary[0] && gmtdefs.verbose) {
		char *type[2] = {"double", "single"};
		fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]);
	}

	if (!y_given) {		/* ensure x-power output */
		output_choice[0] = 'x';
		n_outputs = 1;
	}

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif

	alloc_arrays();

	read_data();

	compute_spectra();

	write_output(output_choice, n_outputs);

	free_space();

	GMT_end (argc, argv);
}

void	alloc_arrays(void)
{
	n_spec = window/2;	/* This means we skip zero frequency; data are detrended  */
	window_2 = 2 * window;		/* This is for complex array stuff  */

	spec = (struct SPEC *) GMT_memory (VNULL, (size_t)n_spec, sizeof(struct SPEC), GMT_program);
	datac = (float *) GMT_memory (VNULL, (size_t)window_2, sizeof(float), GMT_program);
}

void	read_data(void)
{
	int	i, n_alloc, n_expected_fields, n_fields;
	char	buffer[BUFSIZ];
	double	*in;

	if (fp == NULL) {
		fp = GMT_stdin;
#ifdef SET_IO_MODE
		GMT_setmode (0);
#endif
	}
	n_alloc = GMT_CHUNK;
	n_data = 0;
	x = (float *) GMT_memory (VNULL, (size_t)n_alloc, sizeof(float), GMT_program);
	if (y_given) y = (float *)GMT_memory (VNULL, (size_t)n_alloc, sizeof(float), GMT_program);
	n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 2;

	if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (buffer, BUFSIZ, fp);

	while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

		if (GMT_io.status & GMT_IO_MISMATCH) {
			fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected_fields, n_data);
			exit (EXIT_FAILURE);
		}

                x[n_data] = (float)in[0];
		if (y_given) y[n_data] = (float)in[1];
                n_data++;

                if (n_data == n_alloc) {
                        n_alloc += GMT_CHUNK;
			x = (float *) GMT_memory ((void *)x, (size_t)n_alloc, sizeof(float), GMT_program);
			if (y_given) y = (float *)GMT_memory ((void *)y, (size_t)n_alloc, sizeof(float), GMT_program);
                }
        }
        if (fp != GMT_stdin) GMT_fclose(fp);
	x = (float *) GMT_memory ((void *)x, (size_t)n_data, sizeof(float), GMT_program);
	if (y_given) y = (float *)GMT_memory ((void *)y, (size_t)n_data, sizeof(float), GMT_program);
	if (gmtdefs.verbose) fprintf(stderr,"Read %d data points.\n", n_data);
}

void	compute_spectra(void)
{
	int	n_windows, w, i, t_start, t_stop, t, f;
	int	narray;
	float	work = 0.0;
	double	dw, spec_scale, x_varp, y_varp, one_on_nw, co_quad;
	double	xreal, ximag, yreal, yimag, xpower, ypower, co_spec, quad_spec;

	narray = window;

	/* Scale factor for spectral estimates should be 1/4 of amount given in
		Bendat & Piersol eqn 11-102 because I compute 2 * fft in my
		one-sided code below.  However, tests show that I need 1/8 of
		their equation to match variances approximately:  */

	/* This used to read:  spec_scale = 0.5 / (window_2 * dt);  */
	spec_scale = dt / (window_2);

	d_n_windows = (double)n_data / (double)window;

	n_windows = irint (2.0 * d_n_windows) - 1;
	one_on_nw = 1.0 / (double)n_windows;
	dw = (n_windows > 1) ? (double)(n_data - window) / (double)(n_windows - 1) : 1.0;

	for (w = 0; w < n_windows; w++) {
		t_start = (int)floor (0.5 + w * dw);
		t_stop = t_start + window;
		if (y_given) {
			for (t = t_start, i = 0; t < t_stop; t++, i+=2) {
				datac[i] = x[t];
				datac[i+1] = y[t];
			}
		}
		else {
			for (t = t_start, i = 0; t < t_stop; t++, i+=2) {
				datac[i] = x[t];
				datac[i+1] = 0.0;
			}
		}

		detrend_and_hanning();

		GMT_fourt (datac, &narray, 1, -1, 1, &work);

		/* Get one-sided estimates:  */

		x_varp = spec_scale * (datac[0] * datac[0]);
		if (y_given) {
			y_varp = spec_scale * (datac[1] * datac[1]);
			for (i = 0, f = 2; i < n_spec; i++, f+=2) {
				xreal = (i == n_spec - 1) ? datac[f] : datac[f] + datac[window_2 - f];
				ximag = (i == n_spec - 1) ? 0.0 : datac[f+1] - datac[window_2 - f + 1];
				yreal = (i == n_spec - 1) ? datac[f+1] : datac[f+1] + datac[window_2 - f + 1];
				yimag = (i == n_spec - 1) ? 0.0 : datac[window_2 - f] - datac[f];
				xpower = spec_scale * (xreal * xreal + ximag * ximag);
				ypower = spec_scale * (yreal * yreal + yimag * yimag);
				co_spec = spec_scale * (xreal * yreal + ximag * yimag);
				quad_spec = spec_scale * (ximag * yreal - yimag * xreal);

				x_varp += xpower;
				y_varp += ypower;
				spec[i].xpow += xpower;
				spec[i].ypow += ypower;
				/* Temporarily store co-spec in gain:  */
				spec[i].gain += co_spec;
				/* Temporarily store quad-spec in phase:  */
				spec[i].phase += quad_spec;
			}
			x_varp *= (dt/n_spec);
			y_varp *= (dt/n_spec);
		}
		else {
			for (i = 0, f = 2; i < n_spec; i++, f+=2) {
				xreal = datac[f] + datac[window_2 - f];
				ximag = datac[f+1] - datac[window_2 - f + 1];
				xpower = spec_scale * (xreal * xreal + ximag * ximag);
				x_varp += xpower;
				spec[i].xpow += xpower;
			}
			x_varp *= (dt/n_spec);
		}

		if (gmtdefs.verbose) {
			y_pow = (y_given) ? y_variance/y_varp : 0.0;
			fprintf(stderr,"Window %d from %d to %d\n", w, t_start, t_stop);
			sprintf(format, "X var:  %s  X pow:  %s  ratio:  %s  Y var:  %s  Y pow:  %s  ratio:  %s\n",
				gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
			fprintf(stderr, format, x_variance, x_varp, (x_variance/x_varp), y_variance, y_varp, y_pow);
		}
	}
	/* Now we can divide by n_windows for the ensemble average.
		The cross spectral stuff needs to be computed:  */

	if (y_given ) {
		for (i = 0; i < n_spec; i++) {
			spec[i].xpow *= one_on_nw;
			spec[i].ypow *= one_on_nw;
			co_spec = spec[i].gain * one_on_nw;
			quad_spec = spec[i].phase * one_on_nw;
			spec[i].phase = d_atan2(quad_spec, co_spec);
			co_quad = co_spec * co_spec + quad_spec * quad_spec;
			spec[i].coh = co_quad / (spec[i].xpow * spec[i].ypow);
			spec[i].gain = sqrt(co_quad) / spec[i].xpow;
			spec[i].radmit = co_spec / spec[i].xpow;
		}
	}
	else {
		for (i = 0; i < n_spec; i++) {
			spec[i].xpow *= one_on_nw;
		}
	}
}

void	write_output(char *output_choice, int n_outputs)
{
	int	i, j;
	double	delta_f, eps_pow, out[3], *f_or_w;
	char	fname[128];
	FILE	*fpout;
	
	delta_f = 1.0 / (window * dt);
	eps_pow = 1.0 / sqrt(d_n_windows);	/* Multiplicative error bars for power spectra  */

	f_or_w = (double *) GMT_memory (VNULL, (size_t)n_spec, sizeof (double), GMT_program);
	for (i = 0; i < n_spec; i++) f_or_w[i] = (write_wavelength) ? 1.0 / ((i + 1) * delta_f) : (i + 1) * delta_f;

	/* loop through output choices */
	for (j = 0; j < n_outputs; j++) {
	  switch (output_choice[j]) {
	  case 'x':		/* write x power */
	    sprintf(fname, "%s.xpower", namestem);
	    if ( (fpout = GMT_fopen(fname, GMT_io.w_mode)) == NULL) {
	      fprintf(stderr,"%s:  Cannot open w %s\n", GMT_program, fname);
	      exit (EXIT_FAILURE);
	    }
	    if (gmtdefs.verbose) fprintf(stderr,"%s:  Writing %s\n", GMT_program, fname);
	    for (i = 0; i < n_spec; i++) {
	      out[0] = f_or_w[i];
	      out[1] = spec[i].xpow;
	      out[2] = eps_pow * spec[i].xpow;
	      GMT_output (fpout, 3, out);
	    }
	    GMT_fclose(fpout);
	    break;

	  case 'y':		/* Write y power  */
	    sprintf(fname, "%s.ypower", namestem);
	    if ( (fpout = GMT_fopen(fname, GMT_io.w_mode)) == NULL) {
	      fprintf(stderr,"%s:  Cannot open w %s\n", GMT_program, fname);
	      exit (EXIT_FAILURE);
	    }
	    if (gmtdefs.verbose) fprintf(stderr,"%s:  Writing %s\n", GMT_program, fname);
	    for (i = 0; i < n_spec; i++) {
	      out[0] = f_or_w[i];
	      out[1] = spec[i].ypow;
	      out[2] = eps_pow * spec[i].ypow;
	      GMT_output (fpout, 3, out);
	    }
	    GMT_fclose(fpout);
	    break;
	  case 'c':		/* Write Coherent Output power  */
	    sprintf(fname, "%s.cpower", namestem);
	    if ( (fpout = GMT_fopen(fname, GMT_io.w_mode)) == NULL) {
	      fprintf(stderr,"%s:  Cannot open w %s\n", GMT_program, fname);
	      exit (EXIT_FAILURE);
	    }
	    if (gmtdefs.verbose) fprintf(stderr,"%s:  Writing %s\n", GMT_program, fname);
	    for (i = 0; i < n_spec; i++) {
	      out[0] = f_or_w[i];
	      out[1] = spec[i].ypow * spec[i].coh;
	      out[2] = out[1] * eps_pow * sqrt( (2.0 - spec[i].coh) / spec[i].coh);
	      GMT_output (fpout, 3, out);
	    }
	    GMT_fclose(fpout);
	    break;
	  case 'n':		/* Write Noise Output power  */
	    sprintf(fname, "%s.npower", namestem);
	    if ( (fpout = GMT_fopen(fname, GMT_io.w_mode)) == NULL) {
	      fprintf(stderr,"%s:  Cannot open w %s\n", GMT_program, fname);
	      exit (EXIT_FAILURE);
	    }
	    if (gmtdefs.verbose) fprintf(stderr,"%s:  Writing %s\n", GMT_program, fname);
	    for (i = 0; i < n_spec; i++) {
	      out[0] = f_or_w[i];
	      out[1] = spec[i].ypow * (1.0 - spec[i].coh);
	      out[2] = out[1] * eps_pow;
	      GMT_output (fpout, 3, out);
	    }
	    GMT_fclose(fpout);
	    break;
	  case 'g':		/* Write Gain spectrum  */
	    sprintf(fname, "%s.gain", namestem);
	    if ( (fpout = GMT_fopen(fname, GMT_io.w_mode)) == NULL) {
	      fprintf(stderr,"%s:  Cannot open w %s\n", GMT_program, fname);
	      exit (EXIT_FAILURE);
	    }
	    if (gmtdefs.verbose) fprintf(stderr,"%s:  Writing %s\n", GMT_program, fname);
	    for (i = 0; i < n_spec; i++) {
	      out[0] = f_or_w[i];
	      out[1] = spec[i].gain;
	      out[2] = out[1] * eps_pow * sqrt( (1.0 - spec[i].coh) / 2.0 * spec[i].coh);
	      GMT_output (fpout, 3, out);
	    }
	    GMT_fclose(fpout);
	    break;
	  case 'a':		/* Write Real Admittance spectrum  */
	    sprintf(fname, "%s.admit", namestem);
	    if ( (fpout = GMT_fopen(fname, "w")) == NULL) {
	      fprintf(stderr,"%s:  Cannot open w %s\n", GMT_program, fname);
	      exit (EXIT_FAILURE);
	    }
	    if (gmtdefs.verbose) fprintf(stderr,"%s:  Writing %s\n", GMT_program, fname);
	    for (i = 0; i < n_spec; i++) {
	      out[0] = f_or_w[i];
	      out[1] = spec[i].radmit;
	      out[2] = fabs (eps_pow * sqrt( (1.0 - spec[i].coh) / 2.0 * spec[i].coh) * out[1]);
	      GMT_output (fpout, 3, out);
	    }
	    GMT_fclose(fpout);
	    break;
	  case 'p':		/* Write Phase spectrum  */
	    sprintf(fname, "%s.phase", namestem);
	    if ( (fpout = GMT_fopen(fname, GMT_io.w_mode)) == NULL) {
	      fprintf(stderr,"%s:  Cannot open w %s\n", GMT_program, fname);
	      exit (EXIT_FAILURE);
	    }
	    if (gmtdefs.verbose) fprintf(stderr,"%s:  Writing %s\n", GMT_program, fname);
	    for (i = 0; i < n_spec; i++) {
	      out[0] = f_or_w[i];
	      out[1] = spec[i].phase;
	      out[2] = eps_pow * sqrt( (1.0 - spec[i].coh) / 2.0 * spec[i].coh);
	      GMT_output (fpout, 3, out);
	    }
	    GMT_fclose(fpout);
	  case 'o':		/* Write Coherency spectrum  */
	    sprintf(fname, "%s.coh", namestem);
	    if ( (fpout = GMT_fopen(fname, GMT_io.w_mode)) == NULL) {
	      fprintf(stderr,"%s:  Cannot open w %s\n", GMT_program, fname);
	      exit (EXIT_FAILURE);
	    }
	    if (gmtdefs.verbose) fprintf(stderr,"%s:  Writing %s\n", GMT_program, fname);
	    for (i = 0; i < n_spec; i++) {
	      out[0] = f_or_w[i];
	      out[1] = spec[i].coh;
	      out[2] = out[1] * eps_pow * (1.0 - spec[i].coh) * sqrt(2.0 / spec[i].coh);
	      GMT_output (fpout, 3, out);
	    }
	    GMT_fclose(fpout);
	    break;
	  }
	}

	GMT_free ((void *)f_or_w);

}

void	free_space (void)
{
	GMT_free ((void *)spec);
	GMT_free ((void *)datac);
	GMT_free ((void *)x);
	if (y_given) GMT_free ((void *)y);
}

void	detrend_and_hanning(void)
{
	int	i, t;
	double	sumx, sumtx, sumy, sumty, sumt2, x_slope, x_mean, y_slope, y_mean;
	double	t_factor, h_period, h_scale, hc, hw, tt;
	sumx = 0.0;
	sumtx = 0.0;
	sumy = 0.0;
	sumty = 0.0;
	sumt2 = 0.0;
	x_variance = 0.0;
	y_variance = 0.0;
	t_factor = 2.0 / (window - 1);
	h_period = M_PI / (double)window;	/* For Hanning window  */
	h_scale = sqrt(8.0/3.0);		/* For Hanning window  */

	if (y_given) {
		for (i = 0, t = 0; i < window_2; i+=2, t++) {
			tt = t * t_factor - 1.0;
			sumt2 += (tt * tt);
			sumx += datac[i];
			sumtx += (tt * datac[i]);
			sumy += datac[i+1];
			sumty += (tt * datac[i+1]);
		}
	}
	else {
		for (i = 0, t = 0; i < window_2; i+=2, t++) {
			tt = t * t_factor - 1.0;
			sumt2 += (tt * tt);
			sumx += datac[i];
			sumtx += (tt * datac[i]);
		}
	}
	x_slope = sumtx / sumt2;
	x_mean = sumx / window;
	if (y_given) {
		y_slope = sumty / sumt2;
		y_mean = sumy / window;
		for (i = 0, t = 0; i < window_2; i+=2, t++) {
			hc = cos(t * h_period);
			hw = h_scale * (1.0 - hc * hc);
			tt = t * t_factor - 1.0;
			datac[i] -= (float)(x_mean + tt * x_slope);
			datac[i] *= (float)hw;
			x_variance += (datac[i] * datac[i]);
			datac[i+1] -= (float)(y_mean + tt * y_slope);
			datac[i+1] *= (float)hw;
			y_variance += (datac[i+1] * datac[i+1]);
		}
		x_variance /= window;
		y_variance /= window;
	}
	else {
		for (i = 0, t = 0; i < window_2; i+=2, t++) {
			hc = cos(t * h_period);
			hw = h_scale * (1.0 - hc * hc);
			tt = t * t_factor - 1.0;
			datac[i] -= (float)(x_mean + tt * x_slope);
			datac[i] *= (float)hw;
			x_variance += (datac[i] * datac[i]);
		}
		x_variance /= window;
	}
}