home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 2 / goldfish_vol2_cd1.bin / files / misc / sci / gfft / source / okwindow.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-05-25  |  8.7 KB  |  311 lines
  1. /***************************************************************************
  2.  *          Copyright (C) 1994  Charles P. Peterson                  *
  3.  *         4007 Enchanted Sun, San Antonio, Texas 78244-1254             *
  4.  *              Email: Charles_P_Peterson@fcircus.sat.tx.us                *
  5.  *                                                                         *
  6.  *          This is free software with NO WARRANTY.                  *
  7.  *          See gfft.c, or run program itself, for details.              *
  8.  *              Support is available for a fee.                      *
  9.  ***************************************************************************
  10.  *
  11.  * Program:     gfft--General FFT analysis
  12.  * File:        okwindow.c
  13.  * Purpose:     apply window function to a segment of data
  14.  * Author:      Charles Peterson (CPP)
  15.  * History:     30-August-1993 CPP; Created.
  16.  * Comment:     The formulae (and terminology) for the most windows is
  17.  *              from "On the Use of Windows for Harmonic Analysis with the
  18.  *              Discrete Fourier Transform", by Fredric J. Harris,
  19.  *              pp. 172-204 in 'Modern Spectrum Analysis, II', edited
  20.  *              by Stanislav B. Kesler (IEEE Press, New York, 1986).
  21.  *              This classic compendium and comparison of windows was
  22.  *              originally published in 'Proc. IEEE', vol. 66, pp. 51-83,
  23.  *              Jan. 1978.
  24.  *
  25.  *              The formulae for the Parzen, Welch, and Hann[-ing] windows
  26.  *              are taken from 'Numerical Recipes in C', by Press, Flannery,
  27.  *              Teukolsky, and Vetterling (Cambridge University Press,
  28.  *              Cambridge (UK) and New York City, 1988), p. 442.
  29.  *        
  30.  */
  31.  
  32. #include <math.h>
  33.  
  34. #ifdef _AMIGA
  35. #ifdef _M68881
  36. #include <m68881.h>
  37. #endif
  38. #endif
  39.  
  40. #include "gfft.h"
  41. #include "settings.h"
  42.  
  43. #define SQUARE(x) ((x) * (x))
  44.  
  45. void triangle_window (float *vector, ULONG data_count);
  46. void parzen_window (float *vector, ULONG data_count);
  47. void hann_window (float *vector, ULONG data_count);
  48. void hamming_window (float *vector, ULONG data_count);
  49. void welch_window (float *vector, ULONG data_count);
  50. void blackman_harris_74db_window (float *vector, ULONG data_count);
  51. void blackman_harris_92db_window (float *vector, ULONG data_count);
  52. static void blackman_window_n (float *vector, ULONG data_count, 
  53.                    float a0, float a1, float a2, float a3);
  54.  
  55. static float *Window_Vector = NULL;  /* window vector reused if possible */
  56. static unsigned long last_count = 0;
  57. static last_window_type = RECTANGLE_WINDOW;
  58. static double window_gain2 = 1.0;
  59. static double sum_of_window_squares;
  60. static void window_vector__make (ULONG data_count);
  61.  
  62.  
  63. void ok_window (float *indata, ULONG data_count)
  64. {
  65.     ULONG i;
  66. /*
  67.  * Create new window vector if old one...
  68.  *   doesn't exist
  69.  *   is inapplicable
  70.  */
  71.     if (Window_Vector == NULL ||
  72.     last_count != data_count ||
  73.     last_window_type != WindowType)
  74.     {
  75.     window_vector__make (data_count);
  76.     }
  77. /*
  78.  * OK, now apply the window
  79.  */
  80.     if (WindowType == RECTANGLE_WINDOW) return;
  81.     for (i = 0; i < data_count; i++)
  82.     {
  83.     indata[i] *= Window_Vector[i];
  84.     }
  85.     return;
  86. }
  87.  
  88. static void window_vector__make (ULONG data_count)
  89. {
  90.     if (WindowType == RECTANGLE_WINDOW)  /* This is special cased */
  91.     {
  92.     if (Window_Vector) gfree (Window_Vector);
  93.     Window_Vector = NULL;  /* This is VERY important! */
  94.     sum_of_window_squares = data_count;
  95.     window_gain2 = 1.0;
  96.     }
  97.     else
  98.     {
  99.     /*
  100.      * Allocate new memory if
  101.      * none allocated before
  102.      * current allocation is wrong size (in that case, free old memory too)
  103.      */
  104.     if (Window_Vector != NULL && data_count != last_count)
  105.     {
  106.         gfree (Window_Vector);
  107.         Window_Vector = NULL;
  108.         Window_Vector = gmalloc (data_count * sizeof (float),
  109.                      NOTHING_SPECIAL);
  110.     }
  111.     else if (Window_Vector == NULL)
  112.     {
  113.         Window_Vector = gmalloc (data_count * sizeof (float),
  114.                      NOTHING_SPECIAL);
  115.     } /* else, Window_Vector must be the right size already */
  116.  
  117.     sum_of_window_squares = 0.0;
  118.     switch (WindowType)
  119.     {
  120.     case TRIANGLE_WINDOW:
  121.         triangle_window (Window_Vector, data_count);
  122.         break;
  123.     case BLACKMAN_HARRIS_74DB_WINDOW:
  124.         blackman_harris_74db_window (Window_Vector, data_count);
  125.         break;
  126.     case BLACKMAN_HARRIS_92DB_WINDOW:
  127.         blackman_harris_92db_window (Window_Vector, data_count);
  128.         break;
  129.     case HANN_WINDOW:
  130.         hann_window (Window_Vector, data_count);
  131.         break;
  132.     case HAMMING_WINDOW:
  133.         hamming_window (Window_Vector, data_count);
  134.         break;
  135.     case WELCH_WINDOW:
  136.         welch_window (Window_Vector, data_count);
  137.         break;
  138.     case PARZEN_WINDOW:
  139.         parzen_window (Window_Vector, data_count);
  140.         break;
  141.     }
  142.     window_gain2 = sum_of_window_squares / data_count;
  143.     }
  144.     last_count = data_count;
  145.     last_window_type = WindowType;
  146. }
  147.  
  148.  
  149. double ok_window__gain2 (void)
  150. {
  151.     return window_gain2;  /* return hidden value */
  152. }
  153.  
  154.  
  155. void triangle_window (float *vector, ULONG data_count)
  156. {
  157.     ULONG i;
  158.     float half_count = data_count / 2;
  159.  
  160. /*  float half_divisor = (data_count-1.0) / 2.0; */
  161. /*  The above would reach unity at two points nearest center */
  162. /*  The following never reaches unity...unity is at virtual center */
  163. /*  ...which can only be met if window size is odd */
  164. /*  Harris uses the following as well, but see note at end of function */
  165.  
  166.     float half_divisor = (data_count / 2);
  167.  
  168.     for (i = 0; i < half_count; i++)
  169.     {
  170.     vector[i] = i / half_divisor;
  171.     sum_of_window_squares += SQUARE(vector[i]);
  172.     }
  173.     for (i = half_count; i < data_count; i++)
  174.     {
  175.     vector[i] = ((data_count-1) - i) / half_divisor;
  176.     sum_of_window_squares += SQUARE(vector[i]);
  177.     }
  178. /*
  179.  * Based partially on Harris, op. cit., p. 180.
  180.  * But, in the upper part of equation 23b, he appears to be
  181.  * 'off by 1' in the upper range for n, which (I believe) should be:
  182.  * n = 0, 1, ..., (N/2 - 1), as I have done here...
  183.  * He indicates the lower range correctly as beginning with N/2, which
  184.  * is consistent with my interpretation.
  185.  */
  186. }
  187.  
  188. void blackman_harris_74db_window (float *vector, ULONG data_count)
  189. {
  190.     blackman_window_n (vector, data_count,
  191. #ifdef _FFP
  192.                0.40217, 0.49703, 0.09392, 0.00183);
  193. #else
  194.                0.40217F, 0.49703F, 0.09392F, 0.00183F);
  195. #endif
  196. /*
  197.  * Terms from Harris, op. cit., p. 186
  198.  */
  199. }
  200.  
  201. void blackman_harris_92db_window (float *vector, ULONG data_count)
  202. {
  203.     blackman_window_n (vector, data_count,
  204. #ifdef _FFP
  205.                0.35875, 0.48829, 0.14128, 0.001168);
  206. #else
  207.                0.35875F, 0.48829F, 0.14128F, 0.001168F);
  208. #endif
  209. /*
  210.  * Terms from Harris, op. cit., p. 186
  211.  */
  212. }
  213.  
  214. void hamming_window (float *vector, ULONG data_count)
  215. {
  216.     blackman_window_n (vector, data_count,
  217. #ifdef _FFP
  218.                0.54, 0.46, 0.0, 0.0);
  219. #else
  220.                0.54F, 0.46F, 0.0F, 0.0F);
  221. #endif
  222. /*
  223.  * Terms from Harris, op. cit., p. 183
  224.  */
  225. }
  226.  
  227. static void blackman_window_n (float *vector, ULONG data_count, 
  228.                    float a0, float a1, float a2, float a3)
  229. {
  230.     ULONG i;
  231.     float pi_term = 2.0F * PI / data_count;
  232.  
  233.     for (i = 0; i < data_count; i++)
  234.     {
  235.     vector[i] = a0 -
  236.       a1 * cos (pi_term * i) +
  237.         ((a2 != 0.0) ? a2 * cos (pi_term * (2 * i)) : 0.0F) -
  238.           ((a3 != 0.0) ? a3 * cos (pi_term * (3 * i)) : 0.0F);
  239.  
  240.     sum_of_window_squares += SQUARE(vector[i]);
  241.     }
  242. }
  243.  
  244. void parzen_window (float *vector, ULONG data_count)
  245. {
  246.     ULONG i;
  247.  
  248.     for (i = 0; i < data_count; i++)
  249.     {
  250.     vector[i] = 1.0F - (fabs (i - 0.5F * (data_count - 1)) /
  251.                 (0.5F * (data_count + 1)));
  252.  
  253.     sum_of_window_squares += SQUARE(vector[i]);
  254.     }
  255. /*
  256.  * The formula is taken from Press, et. al, op. cit., p. 442.
  257.  * Harris identifies a similarly constructed (but not identical) window
  258.  * as the Riesz window on p. 186.
  259.  */
  260. }
  261.  
  262. void hann_window (float *vector, ULONG data_count)
  263. {
  264.     ULONG i;
  265.  
  266.     for (i = 0; i < data_count; i++)
  267.     {
  268.     vector[i] = 0.5F * (1.0F - cos ( (2.0F * ((float) PI) * i) / 
  269.                        (data_count - 1)));
  270. /*
  271.  * From Press, et. al, op. cit., p 442.
  272.  * In eq. (27b), p. 181, Harris (op. cit.)  seems to be "off by 1" in
  273.  * his divisor, though he gives a footnote identifying the true name,
  274.  * which is used here (not Hanning, which comes from Hann-ing).
  275.  */
  276.  
  277.     sum_of_window_squares += SQUARE(vector[i]);
  278.     }
  279. }
  280.  
  281. void welch_window (float *vector, ULONG data_count)
  282. {
  283.     ULONG i;
  284.     float upper_factor = 0.5F * (data_count - 1);
  285.     float lower_factor = 0.5F * (data_count + 1);
  286.  
  287.     for (i = 0; i < data_count; i++)
  288.     {
  289.     float factor = (i - upper_factor) / lower_factor;
  290. /*
  291.  * From Press, et. al, op. cit., p. 442.
  292.  */
  293.     vector[i] = 1.0F - SQUARE(factor);
  294.  
  295.     sum_of_window_squares += SQUARE(vector[i]);
  296.     }
  297. }
  298.  
  299. /**************** Now, this is special cased away ************************\
  300. void rectangle (float *vector, ULONG data_count)
  301. {
  302.     ULONG i;
  303.  
  304.     for (i = 0; i < data_count; i++)
  305.     {
  306.     vector[i] = 1.0;
  307.     }
  308.     sum_of_window_squares = data_count;
  309. }
  310. \*************************************************************************/
  311.