jDttSP/window.c

   1 /* window.c
   2 This file is part of a program that implements a Software-Defined Radio.
   3
   4 Copyright (C) 2004 by Frank Brickle, AB2KT and Bob McGwier, N4HY
   5 Implemented from code by Bill Schottstaedt of Snd Editor at CCRMA
   6
   7 This program is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU General Public License as published by
   9 the Free Software Foundation; either version 2 of the License, or
  10 (at your option) any later version.
  11
  12 This program is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with this program; if not, write to the Free Software
  19 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  20
  21 The authors can be reached by email at
  22
  23 ab2kt@arrl.net
  24 or
  25 rwmcgwier@comcast.net
  26
  27 or by paper mail at
  28
  29 The DTTS Microwave Society
  30 6 Kathleen Place
  31 Bridgewater, NJ 08807
  32 */
  33
  34 #include <window.h>
  35
  36 /* shamelessly stolen from Bill Schottstaedt's clm.c */
  37 /* made worse in the process, but enough for our purposes here */
  38
  39 //static double
  40 //sqr(double x) { return (x * x); }
  41
  42 /* mostly taken from
  43  *    Fredric J. Harris, "On the Use of Windows for Harmonic Analysis with the
  44  *    Discrete Fourier Transform," Proceedings of the IEEE, Vol. 66, No. 1,
  45  *    January 1978.
  46  *    Albert H. Nuttall, "Some Windows with Very Good Sidelobe Behaviour",
  47  *    IEEE Transactions of Acoustics, Speech, and Signal Processing, Vol. ASSP-29,
  48  *    No. 1, February 1981, pp 84-91
  49  *
  50  * JOS had slightly different numbers for the Blackman-Harris windows.
  51  */
  52
  53 double *
  54 makewindow(Windowtype type, int size, double *window) {
  55   int i, j, midn, midp1, midm1;
  56   double freq, rate, sr1, angle, expn, expsum, cx, two_pi;
  57
  58   midn = size >> 1;
  59   midp1 = (size + 1) / 2;
  60   midm1 = (size - 1) / 2;
  61   two_pi = 8.0 * atan(1.0);
  62   freq = two_pi / (double) size;
  63   rate = 1.0 / (double) midn;
  64   angle = 0.0;
  65   expn = log(2.0) / (double) midn + 1.0;
  66   expsum = 1.0;
  67
  68   switch (type) {
  69   case RECTANGULAR_WINDOW:
  70     for (i = 0; i < size; i++) window[i] = 1.0;
  71     break;
  72   case HANNING_WINDOW:  /* Hann would be more accurate */
  73     for (i = 0, j = size - 1, angle = 0.0; i <= midn; i++, j--, angle += freq)
  74       window[j] = (window[i] = 0.5 - 0.5 * cos(angle));
  75     break;
  76   case WELCH_WINDOW:
  77     for (i = 0, j = size - 1; i <= midn; i++, j--)
  78       window[j] =
  79         (window[i] = 1.0 - sqr((double) (i - midm1) / (double) midp1));
  80     break;
  81   case PARZEN_WINDOW:
  82     for (i = 0, j = size - 1; i <= midn; i++, j--)
  83       window[j] =
  84         (window[i] = 1.0 - fabs((double) (i - midm1) / (double) midp1));
  85     break;
  86   case BARTLETT_WINDOW:
  87     for (i = 0, j = size - 1, angle = 0.0; i <= midn; i++, j--, angle += rate)
  88       window[j] = (window[i] = angle);
  89     break;
  90   case HAMMING_WINDOW:
  91     for (i = 0, j = size - 1, angle = 0.0; i <= midn; i++, j--, angle += freq)
  92       window[j] = (window[i] = 0.54 - 0.46 * cos(angle));
  93     break;
  94   case BLACKMAN2_WINDOW:        /* using Chebyshev polynomial equivalents here */
  95     for (i = 0, j = size - 1, angle = 0.0; i <= midn; i++, j--, angle += freq) {
  96       cx = cos(angle);
  97       window[j] = (window[i] = (.34401 + (cx * (-.49755 + (cx * .15844)))));
  98     }
  99     break;
 100   case BLACKMAN3_WINDOW:
 101     for (i = 0, j = size - 1, angle = 0.0; i <= midn; i++, j--, angle += freq) {
 102       cx = cos(angle);
 103       window[j] =
 104         (window[i] = (.21747 + (cx * (-.45325 + (cx * (.28256 - (cx * .04672)))))));
 105     }
 106     break;
 107   case BLACKMAN4_WINDOW:
 108     for (i = 0, j = size - 1, angle = 0.0; i <= midn; i++, j--, angle += freq) {
 109       cx = cos(angle);
 110       window[j] = (window[i] =
 111                    (.084037 +
 112                     (cx *
 113                      (-.29145 +
 114                       (cx *
 115                        (.375696 + (cx * (-.20762 + (cx * .041194)))))))));
 116     }
 117     break;
 118   case EXPONENTIAL_WINDOW:
 119     for (i = 0, j = size - 1; i <= midn; i++, j--) {
 120       window[j] = (window[i] = expsum - 1.0);
 121       expsum *= expn;
 122     }
 123     break;
 124   case RIEMANN_WINDOW:
 125     sr1 = two_pi / (double) size;
 126     for (i = 0, j = size - 1; i <= midn; i++, j--) {
 127       if (i == midn) window[j] = (window[i] = 1.0);
 128       else {
 129         /* split out because NeXT C compiler can't handle the full expression */
 130         cx = sr1 * (midn - i);
 131         window[i] = sin(cx) / cx;
 132         window[j] = window[i];
 133       }
 134     }
 135     break;
 136    case BLACKMANHARRIS_WINDOW:
 137         {
 138                 double
 139                         a0 = 0.35875,
 140                         a1 = 0.48829,
 141                         a2 = 0.14128,
 142                         a3 = 0.01168;
 143
 144
 145                 for (i = 0; i<size;i++) {
 146                         window[i] = a0 - a1*cos(two_pi*(double)(i+0.5)/(double)size)
 147                                 + a2*cos(2.0*two_pi*(double)(i+0.5)/(double)size)
 148                                 - a3*cos(3.0*two_pi*(double)(i+0.5)/(double)size);
 149                 }
 150         }
 151     break;
 152  default:
 153     return 0;
 154     break;
 155   }
 156
 157   return window;
 158 }