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