1 /***************************************************************************
2 * Copyright (C) 2012 by Pere Ràfols Soler *
3 * sapista2@gmail.com *
4 * *
5 * This program is free software; you can redistribute it and/or modify *
6 * it under the terms of the GNU General Public License as published by *
7 * the Free Software Foundation; either version 2 of the License, or *
8 * (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU General Public License *
16 * along with this program; if not, write to the *
17 * Free Software Foundation, Inc., *
18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19 ***************************************************************************/
20 #ifndef _EQ10Q_FAST_MATH
21 #define _EQ10Q_FAST_MATH
22 #include <math.h>
23 #include <stdlib.h>
24
25 //LUT addr size
26 #define LUT_ADDR_LENGTH 8
27
28 // Compute fast log10 using DSP Guru trick
29 // Observa que akest metode no utilitza per res el signe, de forma que fa el valor absolut de forma intrinseca! t'estalvies el detector de pic!
30 // x: Input value RAW BITS format IEE745
31 // returns: Fast log10 value
fastLog10(int * x_bits,float * Log10LUT)32 static inline float fastLog10(int *x_bits, float *Log10LUT)
33 {
34 // log10(N) = log10(M) + E*log10(2)
35 return (Log10LUT[(*x_bits & 0x7FFFFF) >> (23 - LUT_ADDR_LENGTH)] + (float)((0x0FF & ((*x_bits >> 23))) - 127) * 0.30102999566398119521373889472449f);
36 }
37
38 // Compute fast log using DSP Guru trick
39 // Observa que akest metode no utilitza per res el signe, de forma que fa el valor absolut de forma intrinseca! t'estalvies el detector de pic!
40 // x: Input value RAW BITS format IEE745
41 // returns: Fast log value
fastLog(int * x_bits,float * LogLUT)42 static inline float fastLog(int *x_bits, float *LogLUT)
43 {
44 // log(N) = log(M) + E*log(2)
45 return (LogLUT[(*x_bits & 0x7FFFFF) >> (23 - LUT_ADDR_LENGTH)] + (float)((0x0FF & ((*x_bits >> 23))) - 127) * 0.693147180559945f);
46 }
47
48 // Conversion of binary mantissa 2 decimal
GetBinaryFraction(int x)49 float GetBinaryFraction(int x)
50 {
51 float res = 0;
52 int i;
53 for (i = 22; i >=0; i--)
54 {
55 int b = x >> i;
56 b &= 0x01;
57 res += (float)(b) * powf(2.0f, (float)(i) - 23.0f);
58 }
59
60 return res;
61 }
62
63 // Generate LUT for mantissa log
GenerateLogLUT()64 float* GenerateLogLUT()
65 {
66 int i;
67 float matissa = 1.0f;
68 int size = (int)( powf(2.0f, LUT_ADDR_LENGTH));
69 float *LogLUT = (float *)malloc(sizeof(float)*size);
70
71 for (i = 0; i < size; i++)
72 {
73 int M = (int)i << (23 - LUT_ADDR_LENGTH);
74 matissa = GetBinaryFraction(M) + 1.0f;
75 LogLUT[i] = logf(matissa);
76 }
77 return LogLUT;
78 }
79
80 // Generate LUT for mantissa log10
GenerateLog10LUT()81 float* GenerateLog10LUT()
82 {
83 int i;
84 float matissa = 1.0f;
85 int size = (int)( powf(2.0f, LUT_ADDR_LENGTH));
86 float *Log10LUT = (float *)malloc(sizeof(float)*size);
87
88 for (i = 0; i < size; i++)
89 {
90 int M = (int)i << (23 - LUT_ADDR_LENGTH);
91 matissa = GetBinaryFraction(M) + 1.0f;
92 Log10LUT[i] = log10f(matissa);
93 }
94 return Log10LUT;
95 }
96
97 // Fast dB to Linear conversion using 8 Taylor Series members
98 // Is aprox a 14% faster than Fast_dB2Lin10() but less acurate
99 // Input Range at less 1 dB error: 70dB (from -52 dB to 18 dB)
100 // Input Range at less 0,5 dB error: 63dB (from -51 dB to 12 dB)
Fast_dB2Lin8(float x)101 static inline float Fast_dB2Lin8(float x)
102 {
103 float xa = x + 30.0f; //Error = 0 centrat a -30dB
104
105 //k1 k2 k3 k4 k5
106 //0.031622777f 0.003640707f 0.000209576f 8.04277e-6f 2.3149e-7f
107
108 //k6 k7 k8 k9
109 //5.33025e-9f 1.02278e-10f 1.68217e-12f 2.42083e-14f
110
111 if (x > -60.0f) //Clip on -60dB because the Taylor aprox starts to deviate at these point
112 {
113 return (0.031622777f + xa * (0.003640707f + xa * (0.000209576f + xa * (8.04277e-6f + xa * (2.3149e-7f + xa * (5.33025e-9f + xa * (1.02278e-10f + xa * (1.68217e-12f + xa * 2.42083e-14f))))))));
114 }
115 else
116 {
117 return 0.001f;
118 }
119 }
120
121 // Fast dB to Linear conversion using 10 Taylor Series members
122 // Is aprox a 14% slower than Fast_dB2Lin8() but more acurate
123 // Input Range at less 1 dB error: 89dB (from -67 dB to 22 dB)
124 // Input Range at less 0,5 dB error: 81dB (from -65 dB to 16 dB)
Fast_dB2Lin10(float x)125 static inline float Fast_dB2Lin10(float x)
126 {
127 float xa = x + 40.0f; //Error = 0 centrat a -40dB
128
129 //k1 k2 k3 k4
130 //0.01f 0.001151293f 6.62737e-5f 2.54335e-6f
131
132 //k5 k6 k7 k8
133 //7.32034e-8f 1.68557e-9f 3.23431e-11f 5.31948e-13f
134
135 //k9 k10 k11
136 //7.65535e-15f 9.79283e-17f 1.12744e-18f
137
138 if (x > -67.0f) //Clip on -67dB because the Taylor aprox starts to deviate at these point
139 {
140 return (0.01f + xa * (0.001151293f + xa * (6.62737e-5f + xa * (2.54335e-6f + xa * (7.32034e-8f + xa * (1.68557e-9f + xa * (3.23431e-11f + xa * (5.31948e-13f + xa * (7.65535e-15f + xa * (9.79283e-17f + xa * 1.12744e-18f))))))))));
141 }
142 else
143 {
144 return 0.0004466835921509630484560471330723885330371558666229248046875f;
145 }
146 }
147 #endif