1 /* 2 * Copyright 2019 The Android Open Source Project 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 */ 16 17 #ifndef RESAMPLER_KAISER_WINDOW_H 18 #define RESAMPLER_KAISER_WINDOW_H 19 20 #include <math.h> 21 22 namespace resampler { 23 24 /** 25 * Calculate a Kaiser window centered at 0. 26 */ 27 class KaiserWindow { 28 public: KaiserWindow()29 KaiserWindow() { 30 setStopBandAttenuation(60); 31 } 32 33 /** 34 * @param attenuation typical values range from 30 to 90 dB 35 * @return beta 36 */ setStopBandAttenuation(double attenuation)37 double setStopBandAttenuation(double attenuation) { 38 double beta = 0.0; 39 if (attenuation > 50) { 40 beta = 0.1102 * (attenuation - 8.7); 41 } else if (attenuation >= 21) { 42 double a21 = attenuation - 21; 43 beta = 0.5842 * pow(a21, 0.4) + (0.07886 * a21); 44 } 45 setBeta(beta); 46 return beta; 47 } 48 setBeta(double beta)49 void setBeta(double beta) { 50 mBeta = beta; 51 mInverseBesselBeta = 1.0 / bessel(beta); 52 } 53 54 /** 55 * @param x ranges from -1.0 to +1.0 56 */ operator()57 double operator()(double x) { 58 double x2 = x * x; 59 if (x2 >= 1.0) return 0.0; 60 double w = mBeta * sqrt(1.0 - x2); 61 return bessel(w) * mInverseBesselBeta; 62 } 63 64 // Approximation of a 65 // modified zero order Bessel function of the first kind. 66 // Based on a discussion at: 67 // https://dsp.stackexchange.com/questions/37714/kaiser-window-approximation bessel(double x)68 static double bessel(double x) { 69 double y = cosh(0.970941817426052 * x); 70 y += cosh(0.8854560256532099 * x); 71 y += cosh(0.7485107481711011 * x); 72 y += cosh(0.5680647467311558 * x); 73 y += cosh(0.3546048870425356 * x); 74 y += cosh(0.120536680255323 * x); 75 y *= 2; 76 y += cosh(x); 77 y /= 13; 78 return y; 79 } 80 81 private: 82 double mBeta = 0.0; 83 double mInverseBesselBeta = 1.0; 84 }; 85 86 } // namespace resampler 87 #endif //RESAMPLER_KAISER_WINDOW_H 88