Improvement: approximation of math functions

src/common.h

@@ -24,10 +24,10 @@
 #define __LC3_COMMON_H
 
 #include <lc3.h>
+#include "fastmath.h"
 
 #include <limits.h>
 #include <string.h>
-#include <math.h>
 
 
 /**

src/fastmath.h

@@ -0,0 +1,99 @@
+/******************************************************************************
+ *
+ *  Copyright 2022 Google LLC
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at:
+ *
+ *  http://www.apache.org/licenses/LICENSE-2.0
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ ******************************************************************************/
+
+/**
+ * LC3 - Mathematics function approximation
+ */
+
+#ifndef __LC3_FASTMATH_H
+#define __LC3_FASTMATH_H
+
+#include <math.h>
+
+
+/**
+ * Fast 2^n approximation
+ * x               Operand, range -8 to 8
+ * return          2^x approximation (max relative error ~ 7e-6)
+ */
+static inline float fast_exp2f(float x)
+{
+    float y;
+
+    /* --- Polynomial approx in range -0.5 to 0.5 --- */
+
+    static const float c[] = { 1.27191277e-09, 1.47415221e-07,
+        1.35510312e-05, 9.38375815e-04, 4.33216946e-02 };
+
+    y = (    c[0]) * x;
+    y = (y + c[1]) * x;
+    y = (y + c[2]) * x;
+    y = (y + c[3]) * x;
+    y = (y + c[4]) * x;
+    y = (y + 1.f);
+
+    /* --- Raise to the power of 16  --- */
+
+    y = y*y;
+    y = y*y;
+    y = y*y;
+    y = y*y;
+
+    return y;
+}
+
+/**
+ * Fast log2(x) approximation
+ * x               Operand, greater than 0
+ * return          log2(x) approximation (max absolute error ~ 1e-4)
+ */
+static inline float fast_log2f(float x)
+{
+    float y;
+    int e;
+
+    /* --- Polynomial approx in range 0.5 to 1 --- */
+
+    static const float c[] = {
+        -1.29479677, 5.11769018, -8.42295281, 8.10557963, -3.50567360 };
+
+    x = frexpf(x, &e);
+
+    y = (    c[0]) * x;
+    y = (y + c[1]) * x;
+    y = (y + c[2]) * x;
+    y = (y + c[3]) * x;
+    y = (y + c[4]);
+
+    /* --- Add log2f(2^e) and return --- */
+
+    return e + y;
+}
+
+/**
+ * Fast log10(x) approximation
+ * x               Operand, greater than 0
+ * return          log10(x) approximation (max absolute error ~ 1e-4)
+ */
+static inline float fast_log10f(float x)
+{
+    return log10f(2) * fast_log2f(x);
+}
+
+
+#endif /* __LC3_FASTMATH_H */

src/sns.c

@@ -285,7 +285,7 @@ static void compute_scale_factors(enum lc3_dt dt, enum lc3_srate sr,
     float noise_floor = fmaxf(e_sum * (1e-4f / 64), 0x1p-32f);
 
     for (int i = 0; i < LC3_NUM_BANDS; i++)
-        e[i] = log2f(fmaxf(e[i], noise_floor)) * 0.5f;
+        e[i] = fast_log2f(fmaxf(e[i], noise_floor)) * 0.5f;
 
     /* --- Grouping & scaling --- */
 
@@ -706,7 +706,7 @@ static void spectral_shaping(enum lc3_dt dt, enum lc3_srate sr,
     const int *lim = lc3_band_lim[dt][sr];
 
     for (int i = 0, ib = 0; ib < nb; ib++) {
-        float g_sns = powf(2, -scf[ib]);
+        float g_sns = fast_exp2f(-scf[ib]);
 
         for ( ; i < lim[ib+1]; i++)
             y[i] = x[i] * g_sns;

src/spec.c

@@ -69,7 +69,7 @@ static int estimate_gain(
         x_max = fmaxf(x_max, x3);
 
         float s2 = x0*x0 + x1*x1 + x2*x2 + x3*x3;
-        e[i] = 28.f/20 * 10 * (s2 > 0 ? log10f(s2) : -10);
+        e[i] = 28.f/20 * 10 * (s2 > 0 ? fast_log10f(s2) : -10);
     }
 
     /* --- Determine gain index --- */

tables/fastmath.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+#
+# Copyright 2022 Google LLC
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def fast_exp2(x, p):
+
+    p = p.astype(np.float32)
+    x = x.astype(np.float32)
+
+    y = (((((p[0]*x) + p[1])*x + p[2])*x + p[3])*x + p[4])*x + 1
+
+    return np.power(y.astype(np.float32), 16)
+
+def approx_exp2():
+
+    x = np.arange(-8, 8, step=1e-3)
+
+    p = np.polyfit(x, ((2 ** (x/16)) - 1) / x, 4)
+    y = [ fast_exp2(x[i], p) for i in range(len(x)) ]
+    e = np.abs(y - 2**x) / (2 ** x)
+
+    print('{{ {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e} }}'
+                .format(p[0], p[1], p[2], p[3], p[4]))
+    print('Max relative error: ', np.max(e))
+    print('Max RMS error: ', np.sqrt(np.mean(e ** 2)))
+
+    if False:
+        fig, (ax1, ax2) = plt.subplots(2)
+
+        ax1.plot(x, 2**x, label='Reference')
+        ax1.plot(x, y, label='Approximation')
+        ax1.legend()
+
+        ax2.plot(x, e, label='Relative Error')
+        ax2.legend()
+
+        plt.show()
+
+
+def fast_log2(x, p):
+
+    p = p.astype(np.float32)
+    x = x.astype(np.float32)
+
+    (x, e) = np.frexp(x)
+
+    y = ((((p[0]*x) + p[1])*x + p[2])*x + p[3])*x + p[4]
+
+    return (e ) + y.astype(np.float32)
+
+def approx_log2():
+
+    x = np.logspace(-1, 0, base=2, num=100)
+    p = np.polyfit(x, np.log2(x), 4)
+
+    x = np.logspace(-2, 5, num=10000)
+    y = [ fast_log2(x[i], p) for i in range(len(x)) ]
+    e = np.abs(y - np.log2(x))
+
+    print('{{ {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e}, {:14.8e} }}'
+                .format(p[0], p[1], p[2], p[3], p[4]))
+    print('Max absolute error: ', np.max(e))
+    print('Max RMS error: ', np.sqrt(np.mean(e ** 2)))
+
+    if False:
+        fig, (ax1, ax2) = plt.subplots(2)
+
+        ax1.plot(x, np.log2(x),  label='Reference')
+        ax1.plot(x, y, label='Approximation')
+        ax1.legend()
+
+        ax2.plot(x, e, label = 'Absolute error')
+        ax2.legend()
+
+        plt.show()
+
+
+if __name__ == '__main__':
+
+    print('\n--- Approximation of 2^n ---')
+    approx_exp2()
+
+    print('\n--- Approximation of log2(n) ---')
+    approx_log2()
+
+    print('')

test/sns.py

@@ -477,7 +477,7 @@ def check_analysis(rng, dt, sr):
                 ok = ok and data_c[k] == data[k]
 
             ok = ok and lc3.sns_get_nbits() == analysis.get_nbits()
-            ok = ok and np.amax(np.abs(y - y_c)) < 1e-2
+            ok = ok and np.amax(np.abs(y - y_c)) < 1e-1
 
     return ok
 
@@ -518,7 +518,7 @@ def check_analysis_appendix_c(dt):
     for i in range(len(C.E_B[dt])):
 
         scf = lc3.sns_compute_scale_factors(dt, sr, C.E_B[dt][i], False)
-        ok = ok and np.amax(np.abs(scf - C.SCF[dt][i])) < 1e-5
+        ok = ok and np.amax(np.abs(scf - C.SCF[dt][i])) < 1e-4
 
         (lf, hf) = lc3.sns_resolve_codebooks(scf)
         ok = ok and lf == C.IND_LF[dt][i] and hf == C.IND_HF[dt][i]
@@ -535,7 +535,7 @@ def check_analysis_appendix_c(dt):
         ok = ok and np.amax(np.abs(scf_q - C.SCF_Q[dt][i])) < 1e-5
 
         x = lc3.sns_spectral_shaping(dt, sr, C.SCF_Q[dt][i], False, C.X[dt][i])
-        ok = ok and np.amax(np.abs(1 - x/C.X_S[dt][i])) < 1e-6
+        ok = ok and np.amax(np.abs(1 - x/C.X_S[dt][i])) < 1e-5
 
         (x, data) = lc3.sns_analyze(dt, sr, C.E_B[dt][i], False, C.X[dt][i])
         ok = ok and data['lfcb'] == C.IND_LF[dt][i]
@@ -549,7 +549,7 @@ def check_analysis_appendix_c(dt):
             data['idx_b'] == C.IDX_B[dt][i])
         ok = ok and (C.LS_IND_B[dt][i] is None or
             data['ls_b'] == C.LS_IND_B[dt][i])
-        ok = ok and np.amax(np.abs(1 - x/C.X_S[dt][i])) < 1e-6
+        ok = ok and np.amax(np.abs(1 - x/C.X_S[dt][i])) < 1e-5
 
     return ok