From 168d1648a6f9fac39cb1c53ae2572496adc6b300 Mon Sep 17 00:00:00 2001 From: Seppo Ingalsuo Date: Wed, 27 Jul 2022 12:44:34 +0300 Subject: [PATCH] Math: Add window functions library This patch adds a library for window functions. The window functions are used to mitigate artefacts (from nature of FFT) to spectra when computing FFT for a finite long chunk of audio. The currently provided window functions are Blackman, Hamming, Povey, and rectangular. Signed-off-by: Seppo Ingalsuo --- src/include/sof/math/window.h | 53 +++++++++++++++ src/math/CMakeLists.txt | 4 ++ src/math/Kconfig | 13 ++++ src/math/window.c | 122 ++++++++++++++++++++++++++++++++++ 4 files changed, 192 insertions(+) create mode 100644 src/include/sof/math/window.h create mode 100644 src/math/window.c diff --git a/src/include/sof/math/window.h b/src/include/sof/math/window.h new file mode 100644 index 000000000..bd04317b9 --- /dev/null +++ b/src/include/sof/math/window.h @@ -0,0 +1,53 @@ +/* SPDX-License-Identifier: BSD-3-Clause + * + * Copyright(c) 2022 Intel Corporation. All rights reserved. + * + * Author: Seppo Ingalsuo + */ + +/* Window functions related functions */ + +#ifndef __SOF_MATH_WINDOW_H__ +#define __SOF_MATH_WINDOW_H__ + +#include +#include + +#define WIN_BLACKMAN_A0 Q_CONVERT_FLOAT(7938.0 / 18608.0, 15) /* For "exact" blackman */ + +/** + * \brief Return rectangular window, simply values of one + * \param[in,out] win Output vector with coefficients + * \param[in] length Length of coefficients vector + */ +void win_rectangular_16b(int16_t win[], int length); + +/** + * \brief Calculate Blackman window function, reference + * https://en.wikipedia.org/wiki/Window_function#Blackman_window + * + * \param[in,out] win Output vector with coefficients + * \param[in] length Length of coefficients vector + * \param[in] a0 Parameter for window shape, use e.g. 0.42 as Q1.15 + */ +void win_blackman_16b(int16_t win[], int length, int16_t a0); + +/** + * \brief Calculate Hamming window function, reference + * https://en.wikipedia.org/wiki/Window_function#Hann_and_Hamming_windows + * + * \param[in,out] win Output vector with coefficients + * \param[in] length Length of coefficients vector + */ +void win_hamming_16b(int16_t win[], int length); + +/** + * \brief Calculate Povey window function. It's a window function + * from Pytorch. + * + * \param[in,out] win Output vector with coefficients + * \param[in] length Length of coefficients vector + */ +void win_povey_16b(int16_t win[], int length); + +#endif /* __SOF_MATH_WINDOW_H__ */ diff --git a/src/math/CMakeLists.txt b/src/math/CMakeLists.txt index 54709f9a7..625b1ffa6 100644 --- a/src/math/CMakeLists.txt +++ b/src/math/CMakeLists.txt @@ -45,3 +45,7 @@ endif() if(CONFIG_MATH_IIR_DF2T) add_local_sources(sof iir_df2t_generic.c iir_df2t_hifi3.c iir.c) endif() + +if(CONFIG_MATH_WINDOW) + add_local_sources(sof window.c) +endif() diff --git a/src/math/Kconfig b/src/math/Kconfig index 49e562446..f22f9b2a9 100644 --- a/src/math/Kconfig +++ b/src/math/Kconfig @@ -121,4 +121,17 @@ config MATH_IIR_DF2T Select this to build IIR (Infinite Impulse Response) filter or type 2-transposed library. +config MATH_WINDOW + bool "Window functions library" + default n + select CORDIC_FIXED + select NATURAL_LOGARITHM_FIXED + help + Select this to build a library for window functions. The window functions + are used to mitigate artefacts (from nature of FFT) to spectra when + computing FFT for a finite long chunk of audio. The sample values within + the window length need to be multiplied by the weight of coefficient value. + The currently provided window functions are Blackman, Hamming, Povey, and + rectangular. + endmenu diff --git a/src/math/window.c b/src/math/window.c new file mode 100644 index 000000000..4795112ff --- /dev/null +++ b/src/math/window.c @@ -0,0 +1,122 @@ +// SPDX-License-Identifier: BSD-3-Clause +// +// Copyright(c) 2022 Intel Corporation. All rights reserved. +// +// Author: Seppo Ingalsuo + +#include +#include +#include +#include +#include +#include +#include + +#define WIN_ONE_Q15 INT16_MAX +#define WIN_ONE_Q31 INT32_MAX +#define WIN_05_Q31 Q_CONVERT_FLOAT(0.5, 31) + +#define WIN_TWO_PI_Q28 Q_CONVERT_FLOAT(6.2831853072, 28) + +#define WIN_085_Q31 Q_CONVERT_FLOAT(0.85, 31) +#define WIN_LOG_2POW31_Q26 Q_CONVERT_FLOAT(21.4875625974, 26) + +/* Exact + * #define WIN_HAMMING_A0_Q30 Q_CONVERT_FLOAT(25.0 / 46.0, 30) + * #define WIN_HAMMING_A1_Q30 Q_CONVERT_FLOAT(1 - 25.0 / 46.0, 30) + */ + +/* Common approximations to match e.g. Octave */ +#define WIN_HAMMING_A0_Q30 Q_CONVERT_FLOAT(0.54, 30) +#define WIN_HAMMING_A1_Q30 Q_CONVERT_FLOAT(0.46, 30) + +/** + * \brief Return rectangular window, simply values of one + * \param[in,out] win Output vector with coefficients + * \param[in] length Length of coefficients vector + */ +void win_rectangular_16b(int16_t *win, int length) +{ + int i; + + for (i = 0; i < length; i++) + win[i] = WIN_ONE_Q15; +} + +/** + * \brief Calculate Blackman window function, reference + * https://en.wikipedia.org/wiki/Window_function#Blackman_window + + * \param[in,out] win Output vector with coefficients + * \param[in] length Length of coefficients vector + * \param[in] a0 Parameter for window shape, use e.g. 0.42 as Q1.15 + */ +void win_blackman_16b(int16_t win[], int length, int16_t a0) +{ + const int32_t a1 = Q_CONVERT_FLOAT(0.5, 31); + int32_t inv_length; + int32_t val; + int32_t a; + int16_t alpha; + int32_t a2; + int32_t c1; + int32_t c2; + int n; + + alpha = WIN_ONE_Q15 - 2 * a0; /* Q1.15 */ + a2 = alpha << 15; /* Divided by 2 in Q1.31 */ + a = WIN_TWO_PI_Q28 / (length - 1); /* Q4.28 */ + inv_length = WIN_ONE_Q31 / length; + + for (n = 0; n < length; n++) { + c1 = cos_fixed_32b(a * n); + c2 = cos_fixed_32b(2 * n * Q_MULTSR_32X32((int64_t)a, inv_length, 28, 31, 28)); + val = a0 - Q_MULTSR_32X32((int64_t)a1, c1, 31, 31, 15) + + Q_MULTSR_32X32((int64_t)a2, c2, 31, 31, 15); + win[n] = sat_int16(val); + } +} + +void win_hamming_16b(int16_t win[], int length) +{ + int32_t val; + int32_t a; + int n; + + a = WIN_TWO_PI_Q28 / (length - 1); /* Q4.28 */ + for (n = 0; n < length; n++) { + /* Calculate 0.54 - 0.46 * cos(a * n) */ + val = cos_fixed_32b(a * n); /* Q4.28 -> Q1.31 */ + val = Q_MULTSR_32X32((int64_t)val, WIN_HAMMING_A1_Q30, 31, 30, 30); /* Q1.30 */ + val = WIN_HAMMING_A0_Q30 - val; + + /* Convert to Q1.15 */ + win[n] = sat_int16(Q_SHIFT_RND(val, 30, 15)); + } +} + +void win_povey_16b(int16_t win[], int length) +{ + int32_t cos_an; + int32_t x1; + int32_t x2; + int32_t x3; + int32_t x4; + int32_t a; + int n; + + a = WIN_TWO_PI_Q28 / (length - 1); /* Q4.28 */ + for (n = 0; n < length; n++) { + /* Calculate 0.5 - 0.5 * cos(a * n) */ + cos_an = cos_fixed_32b(a * n); /* Q4.28 -> Q1.31 */ + x1 = WIN_05_Q31 - (cos_an >> 1); /* Q1.31 */ + + /* Calculate x^0.85 as exp(0.85 * log(x)) */ + x2 = (int32_t)(ln_int32((uint32_t)x1) >> 1) - WIN_LOG_2POW31_Q26; + x3 = sat_int32(Q_MULTSR_32X32((int64_t)x2, WIN_085_Q31, 26, 31, 27)); /* Q5.27 */ + x4 = exp_fixed(x3); /* Q5.27 -> Q12.20 */ + + /* Convert to Q1.15 */ + win[n] = sat_int16(Q_SHIFT_RND(x4, 20, 15)); + } +}