/***************************************************************************
* 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;
}
