/***************************************************************************
 *		  Copyright (C) 1994  Charles P. Peterson                  *
 *	     4007 Enchanted Sun, San Antonio, Texas 78244-1254             *
 *              Email: Charles_P_Peterson@fcircus.sat.tx.us                *
 *                                                                         *
 *		  This is free software with NO WARRANTY.                  *
 *	      See gfft.c, or run program itself, for details.              *
 *		      Support is available for a fee.                      *
 ***************************************************************************
 *
 * Program:     gfft--General FFT analysis
 * File:        okwrite.c
 * Purpose:     Write the results of a spectrum analysis
 * Author:      Charles Peterson (CPP)
 * History:     23-August-1993 CPP; Created.
 *              6-Aug-1994 CPP (1.10); Add spectrum file identifier #FFT
 *              28-Aug-94 CPP (1.12); fix LowFrequency for FFP
 *              7-Feb-95 CPP (1.31); Progress Requester
 *              13-Feb-95 CPP (1.40); Header on spectrum files (moved to OK)
 *
 * Comment:     This sure has gotten complicated.  Sorry.
 */

/* #define PARSEVAL_DEBUG */

#define PROGRESS_INCREMENT 128

#include <stdio.h>
#include <math.h>

#ifdef _AMIGA
#ifdef _M68881
#include <m68881.h>
#endif
#endif

#include "gfft.h"
#include "settings.h" /* OkRate, WritePtr, Mean, Amplitude */
#include "errcodes.h"
#include "wbench.h"   /* progress requester */

double LowestFrequencyWritten = 0;
double HighestFrequencyWritten = 0;

void ok_write (BIN_TYPE *bins, long number_bins, long number_segments)
{
    long i;
    double frequency;
    FINAL_TYPE value;
/*    double two_doubles[2]; */
    double nyquist_frequency = OkRate / (2.0L * Interleave);
    double delta_frequency = 0;
    double value_divisor_temp = 1.0 * number_segments;
    FINAL_TYPE multiply = Multiply;
    FINAL_TYPE value_divisor;
    double *mesh;
    FINAL_TYPE smoothing_value_sum = 0.0L;
    double smoothing_frequency_sum = 0.0L;
    FINAL_TYPE save_value;
    double save_frequency;
    double test_frequency;
    long mesh_index = 0;
    long smoothing_count = 0;
    BOOLEAN first_db_of_zero = TRUE;
    double sum_of_bins;
    extern double Ending_Progress;
    double first_write_progress = 100.0 - Ending_Progress;
    double incremental_progress_divisor = number_bins / Ending_Progress;

    if (number_segments == 0) return;

/*
 * reset index used internally in calibration
 * because we're starting from lowest frequency (again, maybe)
 */
    calibration_list__reset (&CalibrationList);

/*
 * If 0 bins (i.e., the DC one only)
 * avoid divide by zero problems
 */
    if (number_bins)
    {
	delta_frequency = nyquist_frequency / number_bins;
	mesh  = ok_mesh (nyquist_frequency, delta_frequency);

	if (Mean)
	{
	/*
         * Mean normalization for power or amplitude.
	 * See Press, et. al, page 439.
	 * To get mean, we'll have to divide by N^2, then by number_segments
         * Note: N is number_bins * 2 if number_bins > 0, thus
	 *  (number_bins * 2)^2 == 4 * number_bins^2
         * 
         * READ THIS: if Parseval option selected (now default)
	 *  we recompute this later using Parseval's theorem.  The
	 *  results are remarkably similar...showing that I WAS doing
	 *  the correct thing here.  The new Parseval method is
	 *  more accurate, but the old method is 2.5% faster overall
	 *  because the squaring and summation in ok_read is not needed.
         */
	    value_divisor_temp = (((4.0 * number_bins) * number_bins) * 
			    number_segments);
	}
	else /* !Mean */
	{

/*          value_divisor_temp = 2.0 * number_bins; */

/*
 * We correct for overlapping bins, by multiplying by extra factor
 * (number_segments * <segment size> / Sample_Frame_Count)
 */

/*	    value_divisor_temp = 2.0 * number_bins * number_segments *
	    2.0 * number_bins / Sample_Frame_Count */
	    
/*
 * which simplifies to the following:
 */

	    value_divisor_temp = 4.0 * number_bins * number_bins *
	      number_segments / Sample_Frame_Count;
	}
    }

    if (WindowType)
/*
 * Special extra correction if non-square window is used
 * NOTE: as currently done, this is an estimate
 * based on window's mean squared average gain...
 * Parseval method works better.
 */
    {
	value_divisor_temp *= ok_window__gain2();
    }

/*
 * This is now where the normalization normally gets done
 * Using Parseval's theorem, i.e.
 *   sum of power in the time domain equals sum of power in freq domain
 */
    if (Parseval)
    {
	sum_of_bins = bins[0];
	if (number_bins > 0)
	{
	    sum_of_bins += bins[number_bins];
	}
	for (i = 1; i <= number_bins-1; i++)
	{
	    sum_of_bins += bins[i] * 2.0;
	}

	if (Mean && Sample_Sum_Of_Squares != 0.0)
	{

#ifdef PARSEVAL_DEBUG
	fprintf (stderr,"Sample MSS: %g, Bin MSS: %g, Old Divisor: %g\n",
		 Sample_Sum_Of_Squares / Sample_Frame_Count,
		 sum_of_bins / value_divisor_temp,
		 value_divisor_temp);
#endif

	    value_divisor_temp = sum_of_bins * Sample_Frame_Count /
	      Sample_Sum_Of_Squares;

#ifdef PARSEVAL_DEBUG
	fprintf (stderr,"Corrected Sample Bin MSS: %g\n",
		 sum_of_bins / value_divisor_temp);
#endif

	}
	else if (!Mean && Sample_Sum_Of_Squares != 0.0)
	{
#ifdef PARSEVAL_DEBUG
	    fprintf (stderr,"Original divisor: %g\n",value_divisor_temp);
#endif

	    value_divisor_temp = sum_of_bins / Sample_Sum_Of_Squares;

#ifdef PARSEVAL_DEBUG
       fprintf (stderr,
		"Post correction: Sample SS: %g, Bin SS: %g, Divisor: %g\n",
		Sample_Sum_Of_Squares,
		sum_of_bins / value_divisor_temp,
		value_divisor_temp);
#endif

	}
	/* else...all's zeros anyway */
    }

    if (PSDensity)
    {
	/*
	 * (Re-)Adjust for PSDensity:
	 * The effect of this is to allow you to splice together curves
	 * computed with different numbers of bins, which is useful in some
	 * applications.
	 *
	 * Ok, we could multiply back in an extra N,
	 * but, I've chosen to divide by the bin size
	 * giving more meaningful units of (.../Hz) or something like that,
	 * I think...
	 */
	value_divisor_temp *= delta_frequency;
    }


    value_divisor = value_divisor_temp;  /* reduce to final precision */

    LowestFrequencyWritten = -1.0;
/*
 * Thanks to bug in FFP handing of DBL_MIN, LOWEST_FREQUENCY is now -1
 * and I have to do this extra stuff
 */
    if (LowFrequency == LOWEST_FREQUENCY)
    {
	i = 1;
    }
    else
    {
	i = 0;
    }
    for (; i <= number_bins; i++)
    {
	if (i % PROGRESS_INCREMENT == 1 && !Time3D)
	{
	    progress_requester__update 
	      ((int) (first_write_progress + 
		      (i / incremental_progress_divisor)));
	}

	frequency = delta_frequency * i;  /* Note 888 below if changing */

	if (frequency < LowFrequency)
	{
	    continue;
	}
	if (frequency > HighFrequency)
	{
	    continue;
	}

	value = bins[i];

	if (i != 0 && i != number_bins)
	{
	    value *= 2.0;    /* Correct for one-sidedness */
	}

	value /= value_divisor;

	if (Pink)
	{
	    if (frequency == 0)
	    {
		continue;
	    }
        /*
         * The following normalizes the spectrum values to about what
         * would be reported by an octave band spectrum analyzer.
         *  ** Warning ** This is a guess, not known to be correct!
         */
	    value *= ((frequency + delta_frequency) / delta_frequency);
	}

	if (Amplitude)
	{
	    value = sqrt (value);
	}

	if (multiply != 1.0)
	{
	    value *= multiply;
	}

	if (mesh && frequency > 0.0L)
	{
	    if (frequency > mesh[mesh_index] && smoothing_count == 0)
	    {
		while (frequency > mesh[++mesh_index]);  /* Catch up */
	    }
	    if (frequency < mesh[mesh_index])
	    {
		smoothing_value_sum += (SquaredSmoothing) ?
		  value * value : value;
		smoothing_frequency_sum += frequency;
		smoothing_count++;
		continue;  /* Nothing to write this time */
	    }
	    else if (frequency == mesh[mesh_index])
	    {
		/* Set up value and frequency for this write */
		smoothing_value_sum += (SquaredSmoothing) ?
		  value * value : value;
		smoothing_frequency_sum += frequency;
		value = smoothing_value_sum / ++smoothing_count;
		if (SquaredSmoothing) value = sqrt (value);
		frequency = smoothing_frequency_sum / smoothing_count;

		/* Set up for next pass */
		smoothing_value_sum = 0.0;
		smoothing_frequency_sum = 0.0;
		smoothing_count = 0;
		mesh_index++;
	    }
	    else if (frequency > mesh[mesh_index])
	    {
		/* Set up value and frequency for this write */
		save_value = value;
		save_frequency = frequency;
		value = smoothing_value_sum / smoothing_count;
		if (SquaredSmoothing) value = sqrt (value);
		frequency = smoothing_frequency_sum / smoothing_count;

		/* Set up for next pass */
		smoothing_value_sum = (SquaredSmoothing) ?
		  save_value * save_value : save_value;
		smoothing_frequency_sum = save_frequency;
		smoothing_count = 1;
		test_frequency = (i+1) * delta_frequency; /*Note 888 above*/
		if (test_frequency > mesh[++mesh_index])
		{   /* Catch up for next pass */
		    while (test_frequency > mesh[++mesh_index]);
		}
	    }
	}		
	if (Db)
	{
	    if (value > 0.0)
	    {
		value = ((Amplitude) ? 20 : 10) * log10 (value);
	    }
	    else
	    {
		if (first_db_of_zero)
		{
		    error_message (DB_OF_ZERO);
		    first_db_of_zero = FALSE;
		}
		continue;  /* Just skip (unless user decides not to) */
	    }
	}

	if (CalibrationList)
	{
	    CATCH_ERROR
	    {
		value = calibration_list__apply (&CalibrationList, 
						 value, frequency);
	    }
	    ON_ERROR
	    {
		if (first_db_of_zero)
		{
		    error_message (DB_OF_ZERO);
		    first_db_of_zero = FALSE;
		}
		continue;  /* Just skip (unless user decides not to) */
	    }
	    END_CATCH_ERROR;
	} /* end if CalibrationList */

	if (Quantization != MIN_QUANTIZATION)
	{
	    long factor = value / Quantization;
	    double value1 = (factor-1) * Quantization;
	    double value2 = factor * Quantization;
	    double value3 = (factor+1) * Quantization;
	    double delta1 = value - value1;
	    double delta2 = fabs (value - value2);
	    double delta3 = value3 - value;
	    if (delta2 <= delta1 && delta2 <= delta3)
	    {
		value = value2;
	    }
	    else if (delta1 <= delta3)
	    {
		value = value1;
	    }
	    else
	    {
		value = value3;
	    }
	}

	if (LowestFrequencyWritten < 0) LowestFrequencyWritten = frequency;
	if (!OutputFormat.binary)
	{
	    if (Time3D)
	    {
#ifdef _FFP
		fprintf (WritePtr, "%-15.8g %-15.8g %-13.7g\n",
			 frequency, OkSegmentTime, value);
#else
		fprintf (WritePtr, "%-19.12g %-19.12g %-13.7g\n", 
			 frequency, OkSegmentTime, value);
#endif
	    }
	    else
	    {
#ifdef _FFP
		fprintf (WritePtr, "%-15.8g %-13.7g\n", frequency, value);
#else
		fprintf (WritePtr, "%-19.12g %-13.7g\n", frequency, value);
#endif
	    }
	}
/*	else  ** This didn't work...need to figure binary formats out **
 *	{
 *	    two_doubles[0] = frequency;
 *	    two_doubles[1] = value;
 *	    fwrite (two_doubles, sizeof(double), (size_t) 2, WritePtr);
 *	}
 */	    
    }
    if (mesh) gfree (mesh);
    if (Time3D) fprintf (WritePtr, "\n");
    HighestFrequencyWritten = frequency;
}