OpenWareLaboratory
KaiserWindow.h
Go to the documentation of this file.
1 /*
2  * Adapted from the Loris C++ Class Library, implementing analysis,
3  * manipulation, and synthesis of digitized sounds using the Reassigned
4  * Bandwidth-Enhanced Additive Sound Model.
5  *
6  * Loris is Copyright (c) 1999-2010 by Kelly Fitz and Lippold Haken
7  *
8  * This program is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY, without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software
20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21  *
22  *
23  * KaiserWindow.C
24  *
25  * Implementation of class Loris::KaiserWindow.
26  *
27  * Kelly Fitz, 14 Dec 1999
28  * loris@cerlsoundgroup.org
29  *
30  * http://www.cerlsoundgroup.org/Loris/
31  *
32  */
33 
34 #include "Window.h"
35 
36 class KaiserWindow : public Window {
37 public:
39  KaiserWindow(float* data, size_t size): Window(data, size){}
40 
41  static KaiserWindow create(float shape, size_t size, bool derivative = false){
42  KaiserWindow win(new float[size], size);
43  if(derivative)
44  buildWindow(win, shape);
45  else
46  buildTimeDerivativeWindow(win, shape);
47  return win;
48  }
49 
50  static void destroy(KaiserWindow obj){
51  Window::destroy(obj);
52  }
53 
54  // ---------------------------------------------------------------------------
55  // buildWindow
56  // ---------------------------------------------------------------------------
68  // α is a non-negative real number that determines the shape of the window.
69  // In the frequency domain, it determines the trade-off between main-lobe width
70  // and side lobe level, which is a central decision in window design.
71  // Sometimes the Kaiser window is parametrized by β, where β = πα.
72  //
73  static void buildWindow(FloatArray win, float shape){
74  // Pre-compute the shared denominator in the Kaiser equation.
75  const float oneOverDenom = 1.0 / zeroethOrderBessel( shape );
76  const unsigned int N = win.getSize() - 1;
77  const float oneOverN = 1.0 / N;
78  float* data = win.getData();
79  for(unsigned int n = 0; n <= N; ++n){
80  const float K = (2.0 * n * oneOverN) - 1.0;
81  const float arg = sqrtf( 1.0 - (K * K) );
82  data[n] = zeroethOrderBessel(shape * arg) * oneOverDenom;
83  // consider using std::cyl_bessel_jf(0, shape * arg) (c++17)
84  }
85  }
86 
87 
88  // ---------------------------------------------------------------------------
89  // createDerivativeWindow
90  // ---------------------------------------------------------------------------
103  //
104  static void buildTimeDerivativeWindow(FloatArray win, float shape){
105  // Pre-compute the common factor that does not depend on n.
106  const unsigned int N = win.getSize() - 1;
107  const float oneOverN = 1.0 / N;
108  const float commonFac = - 2.0 * shape / (N * zeroethOrderBessel(shape) );
109  // w'[0] = w'[N] = 0
110  win[0] = win[N] = 0.0;
111  float* data = win.getData();
112  for(unsigned int n = 1; n < N; ++n){
113  const float K = (2.0 * n * oneOverN) - 1.0;
114  const float arg = sqrtf( 1.0 - (K * K) );
115  data[n] = commonFac * firstOrderBessel(shape * arg) * K / arg;
116  // consider using std::cyl_bessel_jf(1, shape * arg) (c++17)
117  }
118  }
119 
120  // ---------------------------------------------------------------------------
121  // zeroethOrderBessel
122  // ---------------------------------------------------------------------------
123  // Compute the zeroeth order modified Bessel function of the first kind
124  // at x using the series expansion, used to compute the Kasier window
125  // function.
126  //
127  static float zeroethOrderBessel(float x){
128  const float eps = 0.000001;
129 
130  // initialize the series term for m=0 and the result
131  float besselValue = 0;
132  float term = 1;
133  float m = 0;
134 
135  // accumulate terms as long as they are significant
136  while(term > eps * besselValue){
137  besselValue += term;
138 
139  // update the term
140  ++m;
141  term *= (x*x) / (4*m*m);
142  }
143 
144  return besselValue;
145  }
146 
147  // ---------------------------------------------------------------------------
148  // firstOrderBessel
149  // ---------------------------------------------------------------------------
150  // Compute the first order modified Bessel function of the first kind
151  // at x using the series expansion, used to compute the time derivative
152  // of the Kasier window function for computing frequency reassignment.
153  //
154  static float firstOrderBessel(float x){
155  const float eps = 0.000001;
156 
157  // initialize the series term for m=0 and the result
158  float besselValue = 0;
159  float term = .5*x;
160  float m = 0;
161 
162  // accumulate terms as long as they are significant
163  while(term > eps * besselValue)
164  {
165  besselValue += term;
166 
167  // update the term
168  ++m;
169  term *= (x*x) / (4*m*(m+1));
170  }
171 
172  return besselValue;
173  }
174 
175  // ---------------------------------------------------------------------------
176  // computeShape
177  // ---------------------------------------------------------------------------
178  // Compute the Kaiser window shaping parameter from the specified attenuation
179  // of side lobes. This algorithm is given in Kaiser an Schafer,1980 and is
180  // supposed to give better than 0.36% accuracy (Kaiser and Schafer 1980).
181  //
182  static float computeShape(float atten){
183  ASSERT(atten >= 0, "Kaiser window shape must be > 0dB");
184  float alpha;
185  if (atten > 60.0){
186  alpha = 0.12438 * (atten + 6.3);
187  }else if(atten > 13.26){
188  alpha = 0.76609L * ( powf((atten - 13.26), 0.4) ) +
189  0.09834L * (atten - 13.26L);
190  }else{
191  // can't have less than 13dB.
192  alpha = 0.0;
193  }
194  return alpha;
195  }
196 
197  // ---------------------------------------------------------------------------
198  // computeLength
199  // ---------------------------------------------------------------------------
200  // Compute the length (in samples) of the Kaiser window from the desired
201  // (approximate) main lobe width and the control parameter. Of course, since
202  // the window must be an integer number of samples in length, your actual
203  // lobal mileage may vary. This equation appears in Kaiser and Schafer 1980
204  // (on the use of the I0 window class for spectral analysis) as Equation 9.
205  //
206  // The main width of the main lobe must be normalized by the sample rate,
207  // that is, it is a fraction of the sample rate.
208  //
209  static size_t computeLength(float width, float alpha) {
210  // The last 0.5 is cheap rounding.
211  // But I think I don't need cheap rounding because the equation
212  // from Kaiser and Schafer has a +1 that appears to be a cheap
213  // ceiling function.
214  return long(1.0 + (2. * sqrtf((M_PI*M_PI) + (alpha*alpha)) / (M_PI * width)) /* + 0.5 */);
215  }
216 
217 };
#define M_PI
Definition: basicmaths.h:52
This class contains useful methods for manipulating arrays of floats.
Definition: FloatArray.h:12
static void destroy(FloatArray array)
Destroys a FloatArray created with the create() method.
Definition: FloatArray.cpp:472
static void destroy(KaiserWindow obj)
Definition: KaiserWindow.h:50
static void buildWindow(FloatArray win, float shape)
Build a new Kaiser analysis window having the specified shaping parameter.
Definition: KaiserWindow.h:73
static float firstOrderBessel(float x)
Definition: KaiserWindow.h:154
static float computeShape(float atten)
Definition: KaiserWindow.h:182
static void buildTimeDerivativeWindow(FloatArray win, float shape)
Build a new time-derivative Kaiser analysis window having the specified shaping parameter,...
Definition: KaiserWindow.h:104
static size_t computeLength(float width, float alpha)
Definition: KaiserWindow.h:209
KaiserWindow(float *data, size_t size)
Definition: KaiserWindow.h:39
static KaiserWindow create(float shape, size_t size, bool derivative=false)
Definition: KaiserWindow.h:41
static float zeroethOrderBessel(float x)
Definition: KaiserWindow.h:127
size_t getSize() const
Definition: SimpleArray.h:31
T * getData()
Get the data stored in the Array.
Definition: SimpleArray.h:27
Definition: Window.h:20
#define ASSERT(cond, msg)
Definition: message.h:16