/***************************************************************************
 *		  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:        okspctrm.c
 * Purpose:     OK! Do an fft spectrum analysis
 * Author:      Charles Peterson (CPP)
 * History:     24-July-1993 CPP; Created.
 *               6-Feb-95 CPP (1.31); Progress Requester
 *
 * Comment:
 *
 *    As used here, a 'spectrum analysis' is one which starts with real
 *    samples, and ends up with real magnitudes of some sort (possibly with
 *    some sort of 'windowing' applied to the input which may be analyzed
 *    in 'segments' which might be padded and/or overlapped, and possibly
 *    with normalization applied to the output).  The output magnitudes may
 *    either be a 'power' spectrum or an 'amplitude' or 'voltage' spectrum
 *    (the latter being the positive square root of the former.)  *
 *
 *    (Ordinary fft, in contrast, may begin with complex input values and
 *    usually returns complex output values.  Usually the input to an
 *    ordinary fft is not segmented, and doing so may present additional
 *    problems.  Most frequently, ordinary fft is actually part of some
 *    larger problem, such as convolution or spectrum analysis.)
 *
 *    This function has gotten a little out of hand, primarily because of
 *    all the possible combinations of options that are applied here.
 *    Possibly 'Interleave' might be left out next time unless someone
 *    notes a use for it.  (It didn't turn out to be as useful as I
 *    expected.  When curves made with different amounts of interleave are
 *    spliced together, they don't meet, even in 'PSDensity' mode.)  Likewise
 *    for 'Pad,' an option which I've never seen do anything but make
 *    matters much worse.  But somehow, some people seem to think that it's
 *    useful, or at least harmless.  It isn't either, in my experience.
 *    Leaving out those two options would reduce this size of this function
 *    substantially.  
 *
 *    But, don't expect any such streamlining to make a signficant (or
 *    even measureable) difference in performance.  My lprof tests indicate
 *    that this entire beast takes less than 1% of the total time under
 *    typical conditions.  Considering that there is real work going on
 *    here--even if you don't use the more esoteric options--that is quite
 *    reasonable.  (As expected, cfft takes around 70% of the time under
 *    typical conditions, which isn't unreasonable either.)
 *
 *    Anyway, if you want to see my best programming and comment writing,
 *    go to cfft.c, the function which actually does the fft, though it
 *    has gotten harrier since I added an accelerated mode.
 *
 */
#define INITIAL_READ_PROGRESS 5
#define ENDING_WRITE_PROGRESS 45

#define USE_MEMMOVE 
/* Also uses MEMCPY, where applicable */
/* Measured with lprof, only about 0.3% improvement overall. */
/* But, even with 1 bin, didn't make matters worse either. */

#include <time.h>  /* time and difftime */
#include <string.h> /* memmove and memcpy */

#include "gfft.h"
#include "settings.h"
#include "wbench.h"    /* progress requester stuff */

double LoopTimeElapsed = 0.0;

double Percent_Per_Loop = 1.0;
double Initial_Progress = (double) INITIAL_READ_PROGRESS;
double Ending_Progress = (double) ENDING_WRITE_PROGRESS; /* Must not be 0 */

static BIN_TYPE *bins = NULL;
static ULONG number_bins = 0;
static int segment_count = 0;

void do_re_output (void)
{
    if (bins && number_bins && segment_count)
    {
	ok_write (bins, (long) number_bins, (long) segment_count);
    }
    else
    {
	error_message (CANT_RE_OUTPUT);
	RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
    }
}


ULONG ok_spectrum (BOOLEAN do_it_for_real)
{
    static float *indata = NULL;
    static float *overlapdata = NULL;
    static float *interleavedata = NULL;  /* e.g. Aseg Bseg Aseg Bseg ... */
    ULONG actually_read = 0;
    ULONG number_samples = 1;
    ULONG number_interleave_samples = 1;
    ULONG half_samples = 0;
    int more_data = TRUE;
    int error_in_loop = FALSE;
    int overlap = Overlap;     /* May be changed if overlap impossible */
    int last_partial_segment_used = FALSE;
    ULONG loop_count = 0;
    ULONG i;
    ULONG j;
    ULONG ai;
    ULONG aj;
    time_t time_beginning;
    time_t time_ending;
    ULONG total_actually_read;

    number_bins = 0;

    Sample_Sum_Of_Squares = 0.0;
    Sample_Frame_Count = 0;

    total_actually_read = ok_read (NULL, NULL);

    if (NumberBins == MAX_BINS)
    {
	if (!total_actually_read)
	{
	    error_message (NO_DATA_PRESENT);
	    RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
	}
	else
	{
	    ULONG next_power;
	    ULONG total_in_leaf = total_actually_read / Interleave;
	    if (!total_in_leaf)
	    {
		error_message (INSUFFICIENT_DATA);
		RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
	    }
	    next_power = get_pos_power_2 (total_in_leaf);
	    if (Pad || next_power == total_in_leaf)
	    {
		number_bins = next_power / 2; /* pad if required*/
	    }
	    else
	    {
		number_bins = next_power / 4; /* truncate if required */
	    }
	    number_samples = number_bins * 2 * Interleave;
	    bins_d_message (total_actually_read, number_bins);
	}
    }
    else
    {
	number_bins = (ULONG) NumberBins;
	if (number_bins > 0)
	{
	    number_samples = number_bins * 2 * Interleave;
	}
	else
	{
	    number_samples = 1;
	}
    }
/*
 * Now we have number_samples and number_bins
 */
    if (!do_it_for_real) return number_bins;

    segment_count = 0;

    time_beginning = time (NULL);

    half_samples = number_samples / 2;
    if (number_samples < 2)
    {
/*
 * Can't overlap segments of length 1
 */
	overlap = FALSE;
    }
    if (bins)
    {
	gfree (bins);
    }
    if (indata)
    {
	gfree (indata);
    }
    if (overlapdata)
    {
	gfree (overlapdata);
	overlapdata = NULL;
    }
    if (interleavedata)
    {
	gfree (interleavedata);
	interleavedata = NULL;
    }
      
/*
 * Note that there is one extra bin for 0 Hz; bins must be pre-zeroed
 */
    bins = gcalloc ((number_bins + 1), (long) sizeof(BIN_TYPE), 
		    NOTHING_SPECIAL);
    indata = gmalloc (number_samples * sizeof(float), 
		      NOTHING_SPECIAL);
    if (overlap)
    {
	overlapdata = gmalloc (number_samples * sizeof(float), 
			   NOTHING_SPECIAL);
    }
    if (Interleave > 1)
    {
	number_interleave_samples = number_samples / Interleave;
	if (!number_interleave_samples)
	{
	    error_message (INSUFFICIENT_DATA);
	    RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
	}
	interleavedata = gmalloc (number_interleave_samples * sizeof(float),
				 NOTHING_SPECIAL);
    }

/*
 * Now, figure the number of FFT Inner Loops (FILs) required
 *   This is not going to deal with all possible complexities...
 *   since there is no great harm if we are wrong--it's only for the
 *   progress indicator.
 *     In particular, it's not going to deal with Interleave and Pad,
 *     which are fairly useless options not directly supported by the
 *     GUI interface anyway.
 */
    {
	extern long Actual_Time_Segments;
	extern long Inner_Loop_Count;
	long fils_per_segment;
	long est_fils;
	int segments_per_time_segment;
	double computational_percentage;
	int time_segments;

	fils_per_segment = fft_inner_loops (number_bins);

	if (!overlap)
	{
	    segments_per_time_segment = total_actually_read /
	                                 number_samples;
	}
	else if (total_actually_read % number_samples == 0)
	{
	    segments_per_time_segment = 2 * (total_actually_read / 
					 number_samples) - 1;
	}
	else
	{
	    segments_per_time_segment = 2 * (total_actually_read /
					 number_samples);
	}

	if (Time3D)
	{
	    time_segments = Actual_Time_Segments;
	    Ending_Progress = 1;
	    Initial_Progress = 1;
	}
	else
	{
	    Inner_Loop_Count = 0;
	    time_segments = 1;
	    if (segments_per_time_segment * time_segments > 2)
	    {
		Ending_Progress = ENDING_WRITE_PROGRESS * 2 /
		  (segments_per_time_segment * time_segments);
		Initial_Progress = INITIAL_READ_PROGRESS * 2 /
		  (segments_per_time_segment * time_segments);
		if (Ending_Progress < 1) Ending_Progress = 1;
	    }
	    else
	    {
		Ending_Progress = ENDING_WRITE_PROGRESS;
		Initial_Progress = INITIAL_READ_PROGRESS;
	    }
	}
	
	est_fils = fils_per_segment * segments_per_time_segment *
	              time_segments;

	computational_percentage = ((100 - Initial_Progress) -
	                            Ending_Progress);

	if (est_fils)
	{
	    Percent_Per_Loop = computational_percentage / est_fils;
	}
	else
	{
	    Percent_Per_Loop = computational_percentage;
	}
	if (number_bins > 1)
	{
	    progress_requester__update (1);  /* Warm and fuzzy */
	}
	else
	{
	    progress_requester__update (-1);  /* Can't really tell */
	}

    }
/*
 * Now that these preliminaries are out of the way, we can
 * begin looping over input data
 */
    time_beginning = time (NULL);

    while (more_data)
    {
	actually_read = ok_read (indata, number_samples);

	  /* ignored if already > this value */
	progress_requester__update ((int) Initial_Progress);

	if (!actually_read)
	{
	    more_data = FALSE;
	    if (!loop_count)
	    {
		error_message (NO_DATA_PRESENT);
		error_in_loop = TRUE;
	    }
	    break;
	}
	if (actually_read < number_samples)  /* Deal with partial segment */
	{
	    more_data = FALSE;
	    if (overlap && loop_count)
	    {
		for (i = 0, j = actually_read; j < number_samples; 
		     i++, j++)
		{
		    overlapdata[i] = overlapdata[j];
		}
		for (j = 0; j < actually_read; i++, j++)
		{
		    overlapdata[i] = indata[j];
		}
		if (Interleave > 1)    /* Overlap and Interleave */
		{
		    for (ai = 0; ai < Interleave; ai++)
		    {
			float *a_data = interleavedata;
			for (aj = ai; aj < number_samples; aj += Interleave)
			{
			    *a_data++ = overlapdata[aj];
			}
			ok_window (interleavedata, number_interleave_samples);
			ok_rfft (interleavedata, number_interleave_samples);
			ok_sigma (interleavedata, bins, number_bins);
			segment_count++;
		    }
		}
		else  /* Overlap only...no Interleave */
		{
		    ok_window (overlapdata, number_samples);
		    ok_rfft (overlapdata, number_samples);
		    ok_sigma (overlapdata, bins, number_bins);
		    segment_count++;
		}
		overlap = FALSE;  /* overlap buffer now used up */
		last_partial_segment_used = TRUE;
	    }
	    if (Pad)
	    {
		error_message (PADDING_TAILEND);
		for (i = actually_read; i < number_samples; i++)
		{
		    indata[i] = 0.0;
		}
	    }
	    else
	    {
		if (!segment_count)
		{
		    error_message (INSUFFICIENT_DATA);
		    error_in_loop = TRUE;
		}
/*		else if (!last_partial_segment_used)
		{
		    error_message (IGNORING_TAILEND);
		} */
		break;
	    }
	}
/*
 * Here, we've got a full segment to work with (padded or not)
 */
	if (overlap)
	{
	    if (loop_count)  /* Can't process overlap first time around */
	    {
#ifdef USE_MEMMOVE
		memmove (overlapdata, &overlapdata[half_samples], 
			 half_samples * sizeof (float));
		memcpy (&overlapdata[half_samples], indata,
			half_samples * sizeof (float));
#else
		for (i = 0, j = half_samples; i < half_samples; i++, j++)
		{
		    overlapdata[i] = overlapdata[j];
		}
		for (j = 0; j < half_samples; i++, j++)
		{
		    overlapdata[i] = indata[j];
		}
#endif
		if (Interleave > 1)   /* Overlap and Interleave */
		{
		    for (ai = 0; ai < Interleave; ai++)
		    {
			float *a_data = interleavedata;
			for (aj = ai; aj < number_samples; aj += Interleave)
			{
			    *a_data++ = overlapdata[aj];
			}
			ok_window (interleavedata, number_interleave_samples);
			ok_rfft (interleavedata, number_interleave_samples);
			ok_sigma (interleavedata, bins, number_bins);
			segment_count++;
		    }
		}
		else  /* Overlap only...no Interleave */
		{

		    ok_window (overlapdata, number_samples);
		    ok_rfft (overlapdata, number_samples);
		    ok_sigma (overlapdata, bins, number_bins);
		    segment_count++;
		}
	    }
#ifdef USE_MEMMOVE
	    memcpy (overlapdata, indata, number_samples * sizeof (float));
#else
	    for (i = 0; i < number_samples; i++)
	    {
		overlapdata[i] = indata[i];
	    }
#endif
	}
	if (Interleave > 1)  /* Regular (unoverlapped) alternation */
	{
	    for (ai = 0; ai < Interleave; ai++)
	    {
		float *a_data = interleavedata;
		for (aj = ai; aj < number_samples; aj += Interleave)
		{
		    *a_data++ = indata[aj];
		}
		ok_window (interleavedata, number_interleave_samples);
		ok_rfft (interleavedata, number_interleave_samples);
		ok_sigma (interleavedata, bins, number_bins);
		segment_count++;
	    }
	}
	else
	{
	    ok_window (indata, number_samples);
	    ok_rfft (indata, number_samples);
	    ok_sigma (indata, bins, number_bins);
	    segment_count++;
	}
	loop_count++;
    }
    time_ending = time (NULL);
    LoopTimeElapsed = difftime (time_ending, time_beginning);

    if (!error_in_loop)
    {
	loop_time_message (LoopTimeElapsed);
	ok_write (bins, number_bins, segment_count);
    }
    gfree (indata);
    indata = NULL;
/*    gfree (bins); save for ReOut */
/*    bins = NULL; */
    if (overlapdata)
    {
	gfree (overlapdata);
	overlapdata = NULL;
    }
    if (interleavedata)
    {
	gfree (interleavedata);
	interleavedata = NULL;
    }
    if (error_in_loop)
    {
	RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
    }

    return number_bins;
}