mirror of https://github.com/thesofproject/sof.git
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 <seppo.ingalsuo@linux.intel.com>
This commit is contained in:
Seppo Ingalsuo
Liam Girdwood
parent
7288d52ab6
commit
168d1648a6
|
|
@ -0,0 +1,53 @@
|
|||
/* SPDX-License-Identifier: BSD-3-Clause
|
||||
*
|
||||
* Copyright(c) 2022 Intel Corporation. All rights reserved.
|
||||
*
|
||||
* Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
|
||||
*/
|
||||
|
||||
/* Window functions related functions */
|
||||
|
||||
#ifndef __SOF_MATH_WINDOW_H__
|
||||
#define __SOF_MATH_WINDOW_H__
|
||||
|
||||
#include <sof/audio/format.h>
|
||||
#include <stdint.h>
|
||||
|
||||
#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__ */
|
||||
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -0,0 +1,122 @@
|
|||
// SPDX-License-Identifier: BSD-3-Clause
|
||||
//
|
||||
// Copyright(c) 2022 Intel Corporation. All rights reserved.
|
||||
//
|
||||
// Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
|
||||
|
||||
#include <sof/audio/format.h>
|
||||
#include <sof/math/decibels.h>
|
||||
#include <sof/math/log.h>
|
||||
#include <sof/math/trig.h>
|
||||
#include <sof/math/window.h>
|
||||
#include <errno.h>
|
||||
#include <stdint.h>
|
||||
|
||||
#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));
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue