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 #include "ResamplerDefinitions.h"
23 
24 namespace RESAMPLER_OUTER_NAMESPACE::resampler {
25 
26 /**
27  * Calculate a Kaiser window centered at 0.
28  */
29 class KaiserWindow {
30 public:
KaiserWindow()31     KaiserWindow() {
32         setStopBandAttenuation(60);
33     }
34 
35     /**
36      * @param attenuation typical values range from 30 to 90 dB
37      * @return beta
38      */
setStopBandAttenuation(double attenuation)39     double setStopBandAttenuation(double attenuation) {
40         double beta = 0.0;
41         if (attenuation > 50) {
42             beta = 0.1102 * (attenuation - 8.7);
43         } else if (attenuation >= 21) {
44             double a21 = attenuation - 21;
45             beta = 0.5842 * pow(a21, 0.4) + (0.07886 * a21);
46         }
47         setBeta(beta);
48         return beta;
49     }
50 
setBeta(double beta)51     void setBeta(double beta) {
52         mBeta = beta;
53         mInverseBesselBeta = 1.0 / bessel(beta);
54     }
55 
56     /**
57      * @param x ranges from -1.0 to +1.0
58      */
operator()59     double operator()(double x) {
60         double x2 = x * x;
61         if (x2 >= 1.0) return 0.0;
62         double w = mBeta * sqrt(1.0 - x2);
63         return bessel(w) * mInverseBesselBeta;
64     }
65 
66     // Approximation of a
67     // modified zero order Bessel function of the first kind.
68     // Based on a discussion at:
69     // https://dsp.stackexchange.com/questions/37714/kaiser-window-approximation
bessel(double x)70     static double bessel(double x) {
71         double y = cosh(0.970941817426052 * x);
72         y += cosh(0.8854560256532099 * x);
73         y += cosh(0.7485107481711011 * x);
74         y += cosh(0.5680647467311558 * x);
75         y += cosh(0.3546048870425356 * x);
76         y += cosh(0.120536680255323 * x);
77         y *= 2;
78         y += cosh(x);
79         y /= 13;
80         return y;
81     }
82 
83 private:
84     double mBeta = 0.0;
85     double mInverseBesselBeta = 1.0;
86 };
87 
88 } /* namespace RESAMPLER_OUTER_NAMESPACE::resampler */
89 
90 #endif //RESAMPLER_KAISER_WINDOW_H
91