libs/qm-dsp/base/KaiserWindow.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     QM DSP library
   5     Centre for Digital Music, Queen Mary, University of London.
   6
   7     This program is free software; you can redistribute it and/or
   8     modify it under the terms of the GNU General Public License as
   9     published by the Free Software Foundation; either version 2 of the
  10     License, or (at your option) any later version.  See the file
  11     COPYING included with this distribution for more information.
  12 */
  13
  14 #include "KaiserWindow.h"
  15
  16 #include "maths/MathUtilities.h"
  17
  18 KaiserWindow::Parameters
  19 KaiserWindow::parametersForTransitionWidth(double attenuation,
  20                                            double transition)
  21 {
  22     Parameters p;
  23     p.length = 1 + (attenuation > 21.0 ?
  24                     ceil((attenuation - 7.95) / (2.285 * transition)) :
  25                     ceil(5.79 / transition));
  26     p.beta = (attenuation > 50.0 ?
  27               0.1102 * (attenuation - 8.7) :
  28               attenuation > 21.0 ?
  29               0.5842 * pow(attenuation - 21.0, 0.4) + 0.07886 * (attenuation - 21.0) :
  30               0);
  31     return p;
  32 }
  33
  34 static double besselTerm(double x, int i)
  35 {
  36     if (i == 0) {
  37         return 1;
  38     } else {
  39         double f = MathUtilities::factorial(i);
  40         return pow(x/2, i*2) / (f*f);
  41     }
  42 }
  43
  44 static double bessel0(double x)
  45 {
  46     double b = 0.0;
  47     for (int i = 0; i < 20; ++i) {
  48         b += besselTerm(x, i);
  49     }
  50     return b;
  51 }
  52
  53 void
  54 KaiserWindow::init()
  55 {
  56     double denominator = bessel0(m_beta);
  57     bool even = (m_length % 2 == 0);
  58     for (int i = 0; i < (even ? m_length/2 : (m_length+1)/2); ++i) {
  59         double k = double(2*i) / double(m_length-1) - 1.0;
  60         m_window.push_back(bessel0(m_beta * sqrt(1.0 - k*k)) / denominator);
  61     }
  62     for (int i = 0; i < (even ? m_length/2 : (m_length-1)/2); ++i) {
  63         m_window.push_back(m_window[int(m_length/2) - i - 1]);
  64     }
  65 }