blob: 73dbc4178534856df50dec4fcd1a5b4195796b63 [file] [log] [blame]
Phil Burkf0c2a982019-08-17 11:27:37 -07001/*
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
Phil Burk140d7592019-08-21 17:46:59 -070020#include <math.h>
Phil Burkf0c2a982019-08-17 11:27:37 -070021
Phil Burkc123b672019-08-20 10:56:31 -070022namespace resampler {
23
Phil Burkf0c2a982019-08-17 11:27:37 -070024/**
25 * Calculate a Kaiser window centered at 0.
26 */
27class KaiserWindow {
28public:
29 KaiserWindow() {
Phil Burkc123b672019-08-20 10:56:31 -070030 setStopBandAttenuation(60);
Phil Burkf0c2a982019-08-17 11:27:37 -070031 }
32
33 /**
34 * @param attenuation typical values range from 30 to 90 dB
Phil Burkc123b672019-08-20 10:56:31 -070035 * @return beta
Phil Burkf0c2a982019-08-17 11:27:37 -070036 */
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
49 void setBeta(double beta) {
50 mBeta = beta;
Phil Burkc123b672019-08-20 10:56:31 -070051 mInverseBesselBeta = 1.0 / bessel(beta);
Phil Burkf0c2a982019-08-17 11:27:37 -070052 }
53
54 /**
55 * @param x ranges from -1.0 to +1.0
56 */
57 double operator()(double x) {
Phil Burkc123b672019-08-20 10:56:31 -070058 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
Phil Burk140d7592019-08-21 17:46:59 -070064 // Approximation of a
65 // modified zero order Bessel function of the first kind.
Phil Burk6adcdb72019-08-22 17:48:11 -070066 // Based on a discussion at:
Phil Burk140d7592019-08-21 17:46:59 -070067 // https://dsp.stackexchange.com/questions/37714/kaiser-window-approximation
Phil Burkc123b672019-08-20 10:56:31 -070068 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;
Phil Burkf0c2a982019-08-17 11:27:37 -070079 }
80
81private:
82 double mBeta = 0.0;
83 double mInverseBesselBeta = 1.0;
84};
85
Phil Burkc123b672019-08-20 10:56:31 -070086} // namespace resampler
Phil Burkf0c2a982019-08-17 11:27:37 -070087#endif //RESAMPLER_KAISER_WINDOW_H