/***************************************************************************
* 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: okwindow.c
* Purpose: apply window function to a segment of data
* Author: Charles Peterson (CPP)
* History: 30-August-1993 CPP; Created.
* Comment: The formulae (and terminology) for the most windows is
* from "On the Use of Windows for Harmonic Analysis with the
* Discrete Fourier Transform", by Fredric J. Harris,
* pp. 172-204 in 'Modern Spectrum Analysis, II', edited
* by Stanislav B. Kesler (IEEE Press, New York, 1986).
* This classic compendium and comparison of windows was
* originally published in 'Proc. IEEE', vol. 66, pp. 51-83,
* Jan. 1978.
*
* The formulae for the Parzen, Welch, and Hann[-ing] windows
* are taken from 'Numerical Recipes in C', by Press, Flannery,
* Teukolsky, and Vetterling (Cambridge University Press,
* Cambridge (UK) and New York City, 1988), p. 442.
*
*/
#include <math.h>
#ifdef _AMIGA
#ifdef _M68881
#include <m68881.h>
#endif
#endif
#include "gfft.h"
#include "settings.h"
#define SQUARE(x) ((x) * (x))
void triangle_window (float *vector, ULONG data_count);
void parzen_window (float *vector, ULONG data_count);
void hann_window (float *vector, ULONG data_count);
void hamming_window (float *vector, ULONG data_count);
void welch_window (float *vector, ULONG data_count);
void blackman_harris_74db_window (float *vector, ULONG data_count);
void blackman_harris_92db_window (float *vector, ULONG data_count);
static void blackman_window_n (float *vector, ULONG data_count,
float a0, float a1, float a2, float a3);
static float *Window_Vector = NULL; /* window vector reused if possible */
static unsigned long last_count = 0;
static last_window_type = RECTANGLE_WINDOW;
static double window_gain2 = 1.0;
static double sum_of_window_squares;
static void window_vector__make (ULONG data_count);
void ok_window (float *indata, ULONG data_count)
{
ULONG i;
/*
* Create new window vector if old one...
* doesn't exist
* is inapplicable
*/
if (Window_Vector == NULL ||
last_count != data_count ||
last_window_type != WindowType)
{
window_vector__make (data_count);
}
/*
* OK, now apply the window
*/
if (WindowType == RECTANGLE_WINDOW) return;
for (i = 0; i < data_count; i++)
{
indata[i] *= Window_Vector[i];
}
return;
}
static void window_vector__make (ULONG data_count)
{
if (WindowType == RECTANGLE_WINDOW) /* This is special cased */
{
if (Window_Vector) gfree (Window_Vector);
Window_Vector = NULL; /* This is VERY important! */
sum_of_window_squares = data_count;
window_gain2 = 1.0;
}
else
{
/*
* Allocate new memory if
* none allocated before
* current allocation is wrong size (in that case, free old memory too)
*/
if (Window_Vector != NULL && data_count != last_count)
{
gfree (Window_Vector);
Window_Vector = NULL;
Window_Vector = gmalloc (data_count * sizeof (float),
NOTHING_SPECIAL);
}
else if (Window_Vector == NULL)
{
Window_Vector = gmalloc (data_count * sizeof (float),
NOTHING_SPECIAL);
} /* else, Window_Vector must be the right size already */
sum_of_window_squares = 0.0;
switch (WindowType)
{
case TRIANGLE_WINDOW:
triangle_window (Window_Vector, data_count);
break;
case BLACKMAN_HARRIS_74DB_WINDOW:
blackman_harris_74db_window (Window_Vector, data_count);
break;
case BLACKMAN_HARRIS_92DB_WINDOW:
blackman_harris_92db_window (Window_Vector, data_count);
break;
case HANN_WINDOW:
hann_window (Window_Vector, data_count);
break;
case HAMMING_WINDOW:
hamming_window (Window_Vector, data_count);
break;
case WELCH_WINDOW:
welch_window (Window_Vector, data_count);
break;
case PARZEN_WINDOW:
parzen_window (Window_Vector, data_count);
break;
}
window_gain2 = sum_of_window_squares / data_count;
}
last_count = data_count;
last_window_type = WindowType;
}
double ok_window__gain2 (void)
{
return window_gain2; /* return hidden value */
}
void triangle_window (float *vector, ULONG data_count)
{
ULONG i;
float half_count = data_count / 2;
/* float half_divisor = (data_count-1.0) / 2.0; */
/* The above would reach unity at two points nearest center */
/* The following never reaches unity...unity is at virtual center */
/* ...which can only be met if window size is odd */
/* Harris uses the following as well, but see note at end of function */
float half_divisor = (data_count / 2);
for (i = 0; i < half_count; i++)
{
vector[i] = i / half_divisor;
sum_of_window_squares += SQUARE(vector[i]);
}
for (i = half_count; i < data_count; i++)
{
vector[i] = ((data_count-1) - i) / half_divisor;
sum_of_window_squares += SQUARE(vector[i]);
}
/*
* Based partially on Harris, op. cit., p. 180.
* But, in the upper part of equation 23b, he appears to be
* 'off by 1' in the upper range for n, which (I believe) should be:
* n = 0, 1, ..., (N/2 - 1), as I have done here...
* He indicates the lower range correctly as beginning with N/2, which
* is consistent with my interpretation.
*/
}
void blackman_harris_74db_window (float *vector, ULONG data_count)
{
blackman_window_n (vector, data_count,
#ifdef _FFP
0.40217, 0.49703, 0.09392, 0.00183);
#else
0.40217F, 0.49703F, 0.09392F, 0.00183F);
#endif
/*
* Terms from Harris, op. cit., p. 186
*/
}
void blackman_harris_92db_window (float *vector, ULONG data_count)
{
blackman_window_n (vector, data_count,
#ifdef _FFP
0.35875, 0.48829, 0.14128, 0.001168);
#else
0.35875F, 0.48829F, 0.14128F, 0.001168F);
#endif
/*
* Terms from Harris, op. cit., p. 186
*/
}
void hamming_window (float *vector, ULONG data_count)
{
blackman_window_n (vector, data_count,
#ifdef _FFP
0.54, 0.46, 0.0, 0.0);
#else
0.54F, 0.46F, 0.0F, 0.0F);
#endif
/*
* Terms from Harris, op. cit., p. 183
*/
}
static void blackman_window_n (float *vector, ULONG data_count,
float a0, float a1, float a2, float a3)
{
ULONG i;
float pi_term = 2.0F * PI / data_count;
for (i = 0; i < data_count; i++)
{
vector[i] = a0 -
a1 * cos (pi_term * i) +
((a2 != 0.0) ? a2 * cos (pi_term * (2 * i)) : 0.0F) -
((a3 != 0.0) ? a3 * cos (pi_term * (3 * i)) : 0.0F);
sum_of_window_squares += SQUARE(vector[i]);
}
}
void parzen_window (float *vector, ULONG data_count)
{
ULONG i;
for (i = 0; i < data_count; i++)
{
vector[i] = 1.0F - (fabs (i - 0.5F * (data_count - 1)) /
(0.5F * (data_count + 1)));
sum_of_window_squares += SQUARE(vector[i]);
}
/*
* The formula is taken from Press, et. al, op. cit., p. 442.
* Harris identifies a similarly constructed (but not identical) window
* as the Riesz window on p. 186.
*/
}
void hann_window (float *vector, ULONG data_count)
{
ULONG i;
for (i = 0; i < data_count; i++)
{
vector[i] = 0.5F * (1.0F - cos ( (2.0F * ((float) PI) * i) /
(data_count - 1)));
/*
* From Press, et. al, op. cit., p 442.
* In eq. (27b), p. 181, Harris (op. cit.) seems to be "off by 1" in
* his divisor, though he gives a footnote identifying the true name,
* which is used here (not Hanning, which comes from Hann-ing).
*/
sum_of_window_squares += SQUARE(vector[i]);
}
}
void welch_window (float *vector, ULONG data_count)
{
ULONG i;
float upper_factor = 0.5F * (data_count - 1);
float lower_factor = 0.5F * (data_count + 1);
for (i = 0; i < data_count; i++)
{
float factor = (i - upper_factor) / lower_factor;
/*
* From Press, et. al, op. cit., p. 442.
*/
vector[i] = 1.0F - SQUARE(factor);
sum_of_window_squares += SQUARE(vector[i]);
}
}
/**************** Now, this is special cased away ************************\
void rectangle (float *vector, ULONG data_count)
{
ULONG i;
for (i = 0; i < data_count; i++)
{
vector[i] = 1.0;
}
sum_of_window_squares = data_count;
}
\*************************************************************************/
