1 /*
2  * Copyright (C) 2016 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 #include <math.h>
18 #include <audio_utils/limiter.h>
19 
20 #undef USE_ATAN_APPROXIMATION
21 
22 #ifdef USE_ATAN_APPROXIMATION
atan_approximation(float in)23 static inline float atan_approximation(float in)
24 {
25     static const float kPi_2 = 1.5625;  // pi/2 ~= 2.0 * atan_approximation(1.0)
26     if (in <= 1.0) {
27         return in / (1.0 + 0.28 * in * in);
28     } else {
29         const float in_inverse = 1.0 / in;
30         return kPi_2 - in_inverse / (1.0 + 0.28 * in_inverse * in_inverse);
31     }
32 }
33 #endif
34 
limiter(float in)35 float limiter(float in)
36 {
37     static const float crossover = M_SQRT1_2;
38     float in_abs = fabsf(in);
39     if (in_abs <= crossover) {
40         return in;
41     }
42     float out;
43 #ifdef USE_ATAN_APPROXIMATION
44     // atan approximation
45     static const float factor = 4.3125;
46     static const float ratio = 1.0 / factor;
47     out = crossover + atan_approximation((in_abs - crossover)*factor)*ratio;
48     if (out > 1.0) {
49         out = 1.0;
50     }
51 #else
52     // polynomial spline
53     static const float A = 0.3431457505;
54     static const float B = -1.798989873;
55     static const float C = 3.029437252;
56     static const float D = -0.6568542495;
57     // Cubic solution Ax^3 + Bx^2 + Cx + D
58     if (in_abs < M_SQRT2) {
59         // TODO Use Estrin's method P3
60         out = ((A*in_abs + B)*in_abs + C)*in_abs + D;
61     } else {
62         out = 1.0;
63     }
64 #endif
65     if (in < 0) {
66         out = -out;
67     }
68     return out;
69 }
70