libaf/window.c

   1 /*=============================================================================
   2 //
   3 //  This software has been released under the terms of the GNU General Public
   4 //  license. See http://www.gnu.org/copyleft/gpl.html for details.
   5 //
   6 //  Copyright 2001 Anders Johansson ajh@atri.curtin.edu.au
   7 //
   8 //=============================================================================
   9 */
  10
  11 /* Calculates a number of window functions. The following window
  12    functions are currently implemented: Boxcar, Triang, Hanning,
  13    Hamming, Blackman, Flattop and Kaiser. In the function call n is
  14    the number of filter taps and w the buffer in which the filter
  15    coefficients will be stored.
  16 */
  17
  18 #include <math.h>
  19 #include "dsp.h"
  20
  21 /*
  22 // Boxcar
  23 //
  24 // n window length
  25 // w buffer for the window parameters
  26 */
  27 void af_window_boxcar(int n, _ftype_t* w)
  28 {
  29   int i;
  30   // Calculate window coefficients
  31   for (i=0 ; i<n ; i++)
  32     w[i] = 1.0;
  33 }
  34
  35
  36 /*
  37 // Triang a.k.a Bartlett
  38 //
  39 //               |    (N-1)|
  40 //           2 * |k - -----|
  41 //               |      2  |
  42 // w = 1.0 - ---------------
  43 //                    N+1
  44 // n window length
  45 // w buffer for the window parameters
  46 */
  47 void af_window_triang(int n, _ftype_t* w)
  48 {
  49   _ftype_t k1  = (_ftype_t)(n & 1);
  50   _ftype_t k2  = 1/((_ftype_t)n + k1);
  51   int      end = (n + 1) >> 1;
  52   int      i;
  53
  54   // Calculate window coefficients
  55   for (i=0 ; i<end ; i++)
  56     w[i] = w[n-i-1] = (2.0*((_ftype_t)(i+1))-(1.0-k1))*k2;
  57 }
  58
  59
  60 /*
  61 // Hanning
  62 //                   2*pi*k
  63 // w = 0.5 - 0.5*cos(------), where 0 < k <= N
  64 //                    N+1
  65 // n window length
  66 // w buffer for the window parameters
  67 */
  68 void af_window_hanning(int n, _ftype_t* w)
  69 {
  70   int      i;
  71   _ftype_t k = 2*M_PI/((_ftype_t)(n+1)); // 2*pi/(N+1)
  72
  73   // Calculate window coefficients
  74   for (i=0; i<n; i++)
  75     *w++ = 0.5*(1.0 - cos(k*(_ftype_t)(i+1)));
  76 }
  77
  78 /*
  79 // Hamming
  80 //                        2*pi*k
  81 // w(k) = 0.54 - 0.46*cos(------), where 0 <= k < N
  82 //                         N-1
  83 //
  84 // n window length
  85 // w buffer for the window parameters
  86 */
  87 void af_window_hamming(int n,_ftype_t* w)
  88 {
  89   int      i;
  90   _ftype_t k = 2*M_PI/((_ftype_t)(n-1)); // 2*pi/(N-1)
  91
  92   // Calculate window coefficients
  93   for (i=0; i<n; i++)
  94     *w++ = 0.54 - 0.46*cos(k*(_ftype_t)i);
  95 }
  96
  97 /*
  98 // Blackman
  99 //                       2*pi*k             4*pi*k
 100 // w(k) = 0.42 - 0.5*cos(------) + 0.08*cos(------), where 0 <= k < N
 101 //                        N-1                 N-1
 102 //
 103 // n window length
 104 // w buffer for the window parameters
 105 */
 106 void af_window_blackman(int n,_ftype_t* w)
 107 {
 108   int      i;
 109   _ftype_t k1 = 2*M_PI/((_ftype_t)(n-1)); // 2*pi/(N-1)
 110   _ftype_t k2 = 2*k1; // 4*pi/(N-1)
 111
 112   // Calculate window coefficients
 113   for (i=0; i<n; i++)
 114     *w++ = 0.42 - 0.50*cos(k1*(_ftype_t)i) + 0.08*cos(k2*(_ftype_t)i);
 115 }
 116
 117 /*
 118 // Flattop
 119 //                                        2*pi*k                     4*pi*k
 120 // w(k) = 0.2810638602 - 0.5208971735*cos(------) + 0.1980389663*cos(------), where 0 <= k < N
 121 //                                          N-1                        N-1
 122 //
 123 // n window length
 124 // w buffer for the window parameters
 125 */
 126 void af_window_flattop(int n,_ftype_t* w)
 127 {
 128   int      i;
 129   _ftype_t k1 = 2*M_PI/((_ftype_t)(n-1)); // 2*pi/(N-1)
 130   _ftype_t k2 = 2*k1;                   // 4*pi/(N-1)
 131
 132   // Calculate window coefficients
 133   for (i=0; i<n; i++)
 134     *w++ = 0.2810638602 - 0.5208971735*cos(k1*(_ftype_t)i) + 0.1980389663*cos(k2*(_ftype_t)i);
 135 }
 136
 137 /* Computes the 0th order modified Bessel function of the first kind.
 138 // (Needed to compute Kaiser window)
 139 //
 140 // y = sum( (x/(2*n))^2 )
 141 //      n
 142 */
 143 #define BIZ_EPSILON 1E-21 // Max error acceptable
 144
 145 static _ftype_t besselizero(_ftype_t x)
 146 {
 147   _ftype_t temp;
 148   _ftype_t sum   = 1.0;
 149   _ftype_t u     = 1.0;
 150   _ftype_t halfx = x/2.0;
 151   int      n     = 1;
 152
 153   do {
 154     temp = halfx/(_ftype_t)n;
 155     u *=temp * temp;
 156     sum += u;
 157     n++;
 158   } while (u >= BIZ_EPSILON * sum);
 159   return(sum);
 160 }
 161
 162 /*
 163 // Kaiser
 164 //
 165 // n window length
 166 // w buffer for the window parameters
 167 // b beta parameter of Kaiser window, Beta >= 1
 168 //
 169 // Beta trades the rejection of the low pass filter against the
 170 // transition width from passband to stop band.  Larger Beta means a
 171 // slower transition and greater stop band rejection.  See Rabiner and
 172 // Gold (Theory and Application of DSP) under Kaiser windows for more
 173 // about Beta.  The following table from Rabiner and Gold gives some
 174 // feel for the effect of Beta:
 175 //
 176 // All ripples in dB, width of transition band = D*N where N = window
 177 // length
 178 //
 179 // BETA    D       PB RIP   SB RIP
 180 // 2.120   1.50  +-0.27      -30
 181 // 3.384   2.23    0.0864    -40
 182 // 4.538   2.93    0.0274    -50
 183 // 5.658   3.62    0.00868   -60
 184 // 6.764   4.32    0.00275   -70
 185 // 7.865   5.0     0.000868  -80
 186 // 8.960   5.7     0.000275  -90
 187 // 10.056  6.4     0.000087  -100
 188 */
 189 void af_window_kaiser(int n, _ftype_t* w, _ftype_t b)
 190 {
 191   _ftype_t tmp;
 192   _ftype_t k1  = 1.0/besselizero(b);
 193   int      k2  = 1 - (n & 1);
 194   int      end = (n + 1) >> 1;
 195   int      i;
 196
 197   // Calculate window coefficients
 198   for (i=0 ; i<end ; i++){
 199     tmp = (_ftype_t)(2*i + k2) / ((_ftype_t)n - 1.0);
 200     w[end-(1&(!k2))+i] = w[end-1-i] = k1 * besselizero(b*sqrt(1.0 - tmp*tmp));
 201   }
 202 }
 203