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, FLOAT_TYPE* 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, FLOAT_TYPE* w)
  48 {
  49   FLOAT_TYPE k1  = (FLOAT_TYPE)(n & 1);
  50   FLOAT_TYPE k2  = 1/((FLOAT_TYPE)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*((FLOAT_TYPE)(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, FLOAT_TYPE* w)
  69 {
  70   int      i;
  71   FLOAT_TYPE k = 2*M_PI/((FLOAT_TYPE)(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*(FLOAT_TYPE)(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,FLOAT_TYPE* w)
  88 {
  89   int      i;
  90   FLOAT_TYPE k = 2*M_PI/((FLOAT_TYPE)(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*(FLOAT_TYPE)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,FLOAT_TYPE* w)
 107 {
 108   int      i;
 109   FLOAT_TYPE k1 = 2*M_PI/((FLOAT_TYPE)(n-1)); // 2*pi/(N-1)
 110   FLOAT_TYPE 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*(FLOAT_TYPE)i) + 0.08*cos(k2*(FLOAT_TYPE)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,FLOAT_TYPE* w)
 127 {
 128   int      i;
 129   FLOAT_TYPE k1 = 2*M_PI/((FLOAT_TYPE)(n-1)); // 2*pi/(N-1)
 130   FLOAT_TYPE 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*(FLOAT_TYPE)i)
 135                         + 0.1980389663*cos(k2*(FLOAT_TYPE)i);
 136 }
 137
 138 /* Computes the 0th order modified Bessel function of the first kind.
 139 // (Needed to compute Kaiser window)
 140 //
 141 // y = sum( (x/(2*n))^2 )
 142 //      n
 143 */
 144 #define BIZ_EPSILON 1E-21 // Max error acceptable
 145
 146 static FLOAT_TYPE besselizero(FLOAT_TYPE x)
 147 {
 148   FLOAT_TYPE temp;
 149   FLOAT_TYPE sum   = 1.0;
 150   FLOAT_TYPE u     = 1.0;
 151   FLOAT_TYPE halfx = x/2.0;
 152   int      n     = 1;
 153
 154   do {
 155     temp = halfx/(FLOAT_TYPE)n;
 156     u *=temp * temp;
 157     sum += u;
 158     n++;
 159   } while (u >= BIZ_EPSILON * sum);
 160   return(sum);
 161 }
 162
 163 /*
 164 // Kaiser
 165 //
 166 // n window length
 167 // w buffer for the window parameters
 168 // b beta parameter of Kaiser window, Beta >= 1
 169 //
 170 // Beta trades the rejection of the low pass filter against the
 171 // transition width from passband to stop band.  Larger Beta means a
 172 // slower transition and greater stop band rejection.  See Rabiner and
 173 // Gold (Theory and Application of DSP) under Kaiser windows for more
 174 // about Beta.  The following table from Rabiner and Gold gives some
 175 // feel for the effect of Beta:
 176 //
 177 // All ripples in dB, width of transition band = D*N where N = window
 178 // length
 179 //
 180 // BETA    D       PB RIP   SB RIP
 181 // 2.120   1.50  +-0.27      -30
 182 // 3.384   2.23    0.0864    -40
 183 // 4.538   2.93    0.0274    -50
 184 // 5.658   3.62    0.00868   -60
 185 // 6.764   4.32    0.00275   -70
 186 // 7.865   5.0     0.000868  -80
 187 // 8.960   5.7     0.000275  -90
 188 // 10.056  6.4     0.000087  -100
 189 */
 190 void af_window_kaiser(int n, FLOAT_TYPE* w, FLOAT_TYPE b)
 191 {
 192   FLOAT_TYPE tmp;
 193   FLOAT_TYPE k1  = 1.0/besselizero(b);
 194   int      k2  = 1 - (n & 1);
 195   int      end = (n + 1) >> 1;
 196   int      i;
 197
 198   // Calculate window coefficients
 199   for (i=0 ; i<end ; i++){
 200     tmp = (FLOAT_TYPE)(2*i + k2) / ((FLOAT_TYPE)n - 1.0);
 201     w[end-(1&(!k2))+i] = w[end-1-i] = k1 * besselizero(b*sqrt(1.0 - tmp*tmp));
 202   }
 203 }
 204