diff --git a/src/common.h b/src/common.h index daaff5f..751e845 100644 --- a/src/common.h +++ b/src/common.h @@ -24,10 +24,10 @@ #define __LC3_COMMON_H #include +#include "fastmath.h" #include #include -#include /** diff --git a/src/fastmath.h b/src/fastmath.h new file mode 100644 index 0000000..da29024 --- /dev/null +++ b/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 + + +/** + * 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 */ diff --git a/src/sns.c b/src/sns.c index 3f40df5..b06bf60 100644 --- a/src/sns.c +++ b/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; diff --git a/src/spec.c b/src/spec.c index 0f09534..a2c94f8 100644 --- a/src/spec.c +++ b/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 --- */ diff --git a/tables/fastmath.py b/tables/fastmath.py new file mode 100755 index 0000000..9850f6e --- /dev/null +++ b/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('') diff --git a/test/sns.py b/test/sns.py index 01b68e4..897ed79 100644 --- a/test/sns.py +++ b/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