/*************************************************************************** * 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 #ifdef _AMIGA #ifdef _M68881 #include #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; } \*************************************************************************/