diff options
Diffstat (limited to 'src/loris/KaiserWindow.C')
| -rw-r--r-- | src/loris/KaiserWindow.C | 247 |
1 files changed, 247 insertions, 0 deletions
diff --git a/src/loris/KaiserWindow.C b/src/loris/KaiserWindow.C new file mode 100644 index 0000000..c7f3952 --- /dev/null +++ b/src/loris/KaiserWindow.C @@ -0,0 +1,247 @@ +/* + * This is the Loris C++ Class Library, implementing analysis, + * manipulation, and synthesis of digitized sounds using the Reassigned + * Bandwidth-Enhanced Additive Sound Model. + * + * Loris is Copyright (c) 1999-2010 by Kelly Fitz and Lippold Haken + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY, without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * + * KaiserWindow.C + * + * Implementation of class Loris::KaiserWindow. + * + * Kelly Fitz, 14 Dec 1999 + * loris@cerlsoundgroup.org + * + * http://www.cerlsoundgroup.org/Loris/ + * + */ + +#if HAVE_CONFIG_H + #include "config.h" +#endif + +#include "KaiserWindow.h" +#include "LorisExceptions.h" +#include <cmath> + +#if defined(HAVE_M_PI) && (HAVE_M_PI) + const double Pi = M_PI; +#else + const double Pi = 3.14159265358979324; +#endif + +using namespace std; + +// begin namespace +namespace Loris { + +// prototypes for static helpers, defined below +static double zeroethOrderBessel( double x ); +static double firstOrderBessel( double x ); + +// --------------------------------------------------------------------------- +// buildWindow +// --------------------------------------------------------------------------- +//! Build a new Kaiser analysis window having the specified shaping +//! parameter. See Oppenheim and Schafer: "Digital Signal Processing" +//! (1975), p. 452 for further explanation of the Kaiser window. Also, +//! see Kaiser and Schafer, 1980. +//! +//! \param win is the vector that will store the window +//! samples. The number of samples computed will be +//! equal to the length of this vector. Any previous +//! contents will be overwritten. +//! \param shape is the Kaiser shaping parameter, controlling +//! the sidelobe rejection level. +// +void +KaiserWindow::buildWindow( vector< double > & win, double shape ) +{ + // Pre-compute the shared denominator in the Kaiser equation. + const double oneOverDenom = 1.0 / zeroethOrderBessel( shape ); + + const unsigned int N = win.size() - 1; + const double oneOverN = 1.0 / N; + + for ( unsigned int n = 0; n <= N; ++n ) + { + const double K = (2.0 * n * oneOverN) - 1.0; + const double arg = sqrt( 1.0 - (K * K) ); + + win[n] = zeroethOrderBessel( shape * arg ) * oneOverDenom; + } +} + +// --------------------------------------------------------------------------- +// createDerivativeWindow +// --------------------------------------------------------------------------- +//! Build a new time-derivative Kaiser analysis window having the +//! specified shaping parameter, for computing frequency reassignment. +//! The closed form for the time derivative can be obtained from the +//! property of modified Bessel functions that the derivative of the +//! zeroeth order function is equal to the first order function. +//! +//! \param win is the vector that will store the window +//! samples. The number of samples computed will be +//! equal to the length of this vector. Any previous +//! contents will be overwritten. +//! \param shape is the Kaiser shaping parameter, controlling +//! the sidelobe rejection level. +// +void +KaiserWindow::buildTimeDerivativeWindow( vector< double > & win, double shape ) +{ + // Pre-compute the common factor that does not depend on n. + const unsigned int N = win.size() - 1; + const double oneOverN = 1.0 / N; + + const double commonFac = - 2.0 * shape / (N * zeroethOrderBessel( shape ) ); + + // w'[0] = w'[N] = 0 + win[0] = win[N] = 0.0; + + for ( unsigned int n = 1; n < N; ++n ) + { + const double K = (2.0 * n * oneOverN) - 1.0; + const double arg = sqrt( 1.0 - (K * K) ); + + win[n] = commonFac * firstOrderBessel( shape * arg ) * K / arg; + } +} + + +// --------------------------------------------------------------------------- +// zeroethOrderBessel +// --------------------------------------------------------------------------- +// Compute the zeroeth order modified Bessel function of the first kind +// at x using the series expansion, used to compute the Kasier window +// function. +// +static double zeroethOrderBessel( double x ) +{ + const double eps = 0.000001; + + // initialize the series term for m=0 and the result + double besselValue = 0; + double term = 1; + double m = 0; + + // accumulate terms as long as they are significant + while(term > eps * besselValue) + { + besselValue += term; + + // update the term + ++m; + term *= (x*x) / (4*m*m); + } + + return besselValue; +} + +// --------------------------------------------------------------------------- +// firstOrderBessel +// --------------------------------------------------------------------------- +// Compute the first order modified Bessel function of the first kind +// at x using the series expansion, used to compute the time derivative +// of the Kasier window function for computing frequency reassignment. +// +static double firstOrderBessel( double x ) +{ + const double eps = 0.000001; + + // initialize the series term for m=0 and the result + double besselValue = 0; + double term = .5*x; + double m = 0; + + // accumulate terms as long as they are significant + while(term > eps * besselValue) + { + besselValue += term; + + // update the term + ++m; + term *= (x*x) / (4*m*(m+1)); + } + + return besselValue; +} + +// --------------------------------------------------------------------------- +// computeShape +// --------------------------------------------------------------------------- +// Compute the Kaiser window shaping parameter from the specified attenuation +// of side lobes. This algorithm is given in Kaiser an Schafer,1980 and is +// supposed to give better than 0.36% accuracy (Kaiser and Schafer 1980). +// +double +KaiserWindow::computeShape( double atten ) +{ + if ( atten < 0. ) + { + Throw( InvalidArgument, + "Kaiser window shape must be computed from positive (> 0dB)" + " sidelobe attenuation. (received attenuation < 0)" ); + } + + double alpha; + + if ( atten > 60.0 ) + { + alpha = 0.12438 * (atten + 6.3); + } + else if ( atten > 13.26 ) + { + alpha = 0.76609L * ( pow((atten - 13.26), 0.4) ) + + 0.09834L * (atten - 13.26L); + } + else + { + // can't have less than 13dB. + alpha = 0.0; + } + + return alpha; +} +// --------------------------------------------------------------------------- +// computeLength +// --------------------------------------------------------------------------- +// Compute the length (in samples) of the Kaiser window from the desired +// (approximate) main lobe width and the control parameter. Of course, since +// the window must be an integer number of samples in length, your actual +// lobal mileage may vary. This equation appears in Kaiser and Schafer 1980 +// (on the use of the I0 window class for spectral analysis) as Equation 9. +// +// The main width of the main lobe must be normalized by the sample rate, +// that is, it is a fraction of the sample rate. +// +unsigned long +KaiserWindow::computeLength( double width, double alpha ) +{ + //double alpha = computeShape( atten ); + + // The last 0.5 is cheap rounding. + // But I think I don't need cheap rounding because the equation + // from Kaiser and Schafer has a +1 that appears to be a cheap + // ceiling function. + return long(1.0 + (2. * sqrt((Pi*Pi) + (alpha*alpha)) / (Pi * width)) /* + 0.5 */); +} + +} // end of namespace Loris + |