libaf/window.c

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