From 105e514302cba1a20f300657cde060c25c013ecc Mon Sep 17 00:00:00 2001 From: amishutin Date: Thu, 11 Jan 2024 14:43:22 +0300 Subject: [PATCH] added signal module --- modules/README.md | 2 + modules/signal/CMakeLists.txt | 2 + modules/signal/README.md | 4 + modules/signal/include/opencv2/signal.hpp | 17 + .../opencv2/signal/signal_resample.hpp | 32 ++ modules/signal/perf/perf_main.cpp | 6 + modules/signal/perf/perf_precomp.hpp | 15 + modules/signal/perf/perf_resample.cpp | 37 ++ modules/signal/src/precomp.hpp | 10 + modules/signal/src/signal_resample.cpp | 375 ++++++++++++++++++ modules/signal/test/test_main.cpp | 6 + modules/signal/test/test_precomp.hpp | 10 + modules/signal/test/test_signal.cpp | 180 +++++++++ 13 files changed, 696 insertions(+) create mode 100644 modules/signal/CMakeLists.txt create mode 100644 modules/signal/README.md create mode 100644 modules/signal/include/opencv2/signal.hpp create mode 100644 modules/signal/include/opencv2/signal/signal_resample.hpp create mode 100644 modules/signal/perf/perf_main.cpp create mode 100644 modules/signal/perf/perf_precomp.hpp create mode 100644 modules/signal/perf/perf_resample.cpp create mode 100644 modules/signal/src/precomp.hpp create mode 100644 modules/signal/src/signal_resample.cpp create mode 100644 modules/signal/test/test_main.cpp create mode 100644 modules/signal/test/test_precomp.hpp create mode 100644 modules/signal/test/test_signal.cpp diff --git a/modules/README.md b/modules/README.md index cd8ea0fbe..413523f7d 100644 --- a/modules/README.md +++ b/modules/README.md @@ -72,6 +72,8 @@ $ cmake -D OPENCV_EXTRA_MODULES_PATH=/modules -D BUILD_opencv_ + +namespace cv { +namespace signal { + +//! @addtogroup signal +//! @{ + +/** @brief Signal resampling + * + * @param[in] inputSignal Array with input signal. + * @param[out] outSignal Array with output signal + * @param[in] inFreq Input signal frequency. + * @param[in] outFreq Output signal frequency. + * Signal resampling implemented a cubic interpolation function and a filtering function based on Kaiser window and Bessel function, used to construct a FIR filter. + * Result is similar to `scipy.signal.resample`. + +Detail: https://en.wikipedia.org/wiki/Sample-rate_conversion +*/ +CV_EXPORTS_W void resampleSignal(InputArray inputSignal, OutputArray outSignal, const int inFreq, const int outFreq); + +//! @} + +} +} +#endif diff --git a/modules/signal/perf/perf_main.cpp b/modules/signal/perf/perf_main.cpp new file mode 100644 index 000000000..442e2e430 --- /dev/null +++ b/modules/signal/perf/perf_main.cpp @@ -0,0 +1,6 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. +#include "perf_precomp.hpp" + +CV_PERF_TEST_MAIN(signal) diff --git a/modules/signal/perf/perf_precomp.hpp b/modules/signal/perf/perf_precomp.hpp new file mode 100644 index 000000000..2dc91dc21 --- /dev/null +++ b/modules/signal/perf/perf_precomp.hpp @@ -0,0 +1,15 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. +#ifndef __OPENCV_PERF_PRECOMP_HPP__ +#define __OPENCV_PERF_PRECOMP_HPP__ + +#include "opencv2/ts.hpp" +#include "opencv2/signal.hpp" + +namespace opencv_test { +using namespace perf; +using namespace cv::signal; +} + +#endif diff --git a/modules/signal/perf/perf_resample.cpp b/modules/signal/perf/perf_resample.cpp new file mode 100644 index 000000000..b79b7c420 --- /dev/null +++ b/modules/signal/perf/perf_resample.cpp @@ -0,0 +1,37 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. + +#include "perf_precomp.hpp" + +using namespace std; +using namespace cv; +using namespace perf; + +namespace opencv_test { namespace { + +typedef TestBaseWithParam< tuple > TestResampleFunc; + +PERF_TEST_P( TestResampleFunc, resample_sin_signal, + testing::Combine( + testing::Values(1234U, 12345U, 123456U, 1234567U, 12345678U), + testing::Values(16000U, 32000U, 44100U, 48000U), + testing::Values(48000U, 44100U, 32000U, 16000U)) +) +{ + uint32_t sample_signal_size = GET_PARAM(0); + uint32_t inFreq = GET_PARAM(1); + uint32_t outFreq = GET_PARAM(2); + + Mat1f sample_signal(Size(sample_signal_size,1U)); + Mat1f outSignal(Size(1U, 1U)); + for (uint32_t i = 0U; i < (uint32_t)sample_signal.cols; ++i) + { + sample_signal.at(0, i) = sinf(float(i)); + } + declare.in(sample_signal).out(outSignal); + TEST_CYCLE() resampleSignal(sample_signal, outSignal, inFreq, outFreq); + SANITY_CHECK_NOTHING(); +} + +}} diff --git a/modules/signal/src/precomp.hpp b/modules/signal/src/precomp.hpp new file mode 100644 index 000000000..6c2978318 --- /dev/null +++ b/modules/signal/src/precomp.hpp @@ -0,0 +1,10 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html + +#ifndef __OPENCV_SIGNAL_PRECOMP__ +#define __OPENCV_SIGNAL_PRECOMP__ + +#include + +#endif diff --git a/modules/signal/src/signal_resample.cpp b/modules/signal/src/signal_resample.cpp new file mode 100644 index 000000000..03fc5cad3 --- /dev/null +++ b/modules/signal/src/signal_resample.cpp @@ -0,0 +1,375 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html + + +#include "precomp.hpp" + +#include +#include +#include +#include + +#include +#include +#include + +namespace cv { +namespace signal { + +#if (CV_SIMD || CV_SIMD_SCALABLE) +#define v_float32_width (uint32_t)VTraits::vlanes() +const uint32_t v_float32_max_width = (uint32_t)VTraits::max_nlanes; +#endif + +// Modified Bessel function 1st kind 0th order +static float Bessel(float x) +{ + int k = 12; // approximation parameter + float defmul = x * x * 0.25f; + float mul = defmul; + float acc = 0.f; + for(int i = 0 ; i < k; ++i) + { + mul = powf(defmul, static_cast(i)); + mul = mul / powf(tgammaf(static_cast(i + 1)), 2.f); // tgamma(i+1) equals i! + acc +=mul; + } + return acc; +} + +static void init_filter(float beta, int ntabs, float* tabs) +{ + float fc = 0.25f; + // build sinc filter + for (int i = 0; i < ntabs; ++i) + { + tabs[i] = 2 * fc * (i - (ntabs - 1) / 2); + } + std::vector tmparr(ntabs); + for (int i = 0 ; i < ntabs; ++i) + { + if (tabs[i] == 0.f) + { + tmparr[i] = 1.f; + continue; + } + tmparr[i] = (float)(CV_PI * tabs[i]); + } + float mult = 2.f / (float)(ntabs - 1); + // multiply by Kaiser window + for (int i = 0; i < ntabs; ++i) + { + tabs[i] = std::sin(tmparr[i]) / tmparr[i]; + tabs[i] *= Bessel(beta * sqrtf((float)1 - powf((i * mult - 1), 2))) / Bessel(beta); + } + float sum = 0.f; + for (int i = 0 ; i < ntabs; ++i) + { + sum += tabs[i]; + } + sum = 1.f/sum; + // normalize tabs to get unity gain + for (int i = 0; i < ntabs; ++i) + { + tabs[i] *= sum; + } +} + +/////////////// cubic Hermite spline (tail of execIntrinLoop or scalar version) /////////////// +static float scal_cubicHermite(float A, float B, float C, float D, float t) +{ + float a = (-A + (3.0f * B) - (3.0f * C) + D) * 0.5f; + float b = A + C + C - (5.0f * B + D) * 0.5f; + float c = (-A + C) * 0.5f; + return a * t * t * t + b * t * t + c * t + B; +} + +/////////////// cubic Hermite spline (OpenCV's Universal Intrinsics) /////////////// +#if (CV_SIMD || CV_SIMD_SCALABLE) +static inline v_float32 simd_cubicHermite(const v_float32 &v_A, const v_float32 &v_B, const v_float32 &v_C, + const v_float32 &v_D, const v_float32 &v_t) +{ + v_float32 v_zero = vx_setzero_f32(); + v_float32 v_three= vx_setall_f32(3.0f); + v_float32 v_half = vx_setall_f32(0.5f); + v_float32 v_five = vx_setall_f32(5.0f); + + v_float32 v_inv_A = v_sub(v_zero, v_A); + + v_float32 v_a = v_mul(v_sub(v_fma(v_three, v_B, v_add(v_inv_A, v_D)), v_mul(v_three, v_C)), v_half); + v_float32 v_b = v_sub(v_add(v_A, v_C, v_C), v_mul(v_fma(v_five, v_B, v_D), v_half)); + v_float32 v_c = v_mul(v_add(v_inv_A, v_C), v_half); + + return v_add(v_mul(v_a, v_t, v_t, v_t), v_mul(v_b, v_t, v_t), v_fma(v_c, v_t, v_B)); +} +#endif + +static void cubicInterpolate(const Mat1f &src, uint32_t dstlen, Mat1f &dst, uint32_t srclen) +{ + Mat1f tmp(Size(srclen + 3U, 1U)); + tmp.at(0) = src.at(0); + +#if (CV_SIMD || CV_SIMD_SCALABLE) + v_float32 v_reg = vx_setall_f32(src.at(srclen - 1U)); + vx_store(tmp.ptr(0) + (srclen - 1U), v_reg); +#else // scalar version + tmp.at(srclen + 1U) = src.at(srclen - 1U); + tmp.at(srclen + 2U) = src.at(srclen - 1U); +#endif + + uint32_t i = 0U; + +#if (CV_SIMD || CV_SIMD_SCALABLE) + uint32_t len_sub_vfloatStep = (uint32_t)std::max((int64_t)srclen - (int64_t)v_float32_width, (int64_t)0); + for (; i < len_sub_vfloatStep; i+= v_float32_width) + { + v_float32 v_copy = vx_load(src.ptr(0) + i); + vx_store(tmp.ptr(0) + (i + 1U), v_copy); + } +#endif + + // if the tail exists or scalar version + for (; i < srclen; ++i) + { + tmp.at(i + 1U) = src.at(i); + } + + i = 0U; + +#if (CV_SIMD || CV_SIMD_SCALABLE) + int ptr_x_int[v_float32_max_width]; + uint32_t j; + + v_float32 v_dstlen_sub_1 = vx_setall_f32((float)(dstlen - 1U)); + v_float32 v_one = vx_setall_f32(1.0f); + v_float32 v_x_start = v_div(v_one, v_dstlen_sub_1); + v_float32 v_u = vx_setall_f32((float)srclen); + v_float32 v_half = vx_setall_f32(0.5f); + + len_sub_vfloatStep = (uint32_t)std::max((int64_t)dstlen - (int64_t)v_float32_width, (int64_t)0); + for (; i < v_float32_width; ++i) + { + ptr_x_int[i] = (int)i; + } + + float ptr_for_cubicHermite[v_float32_max_width]; + v_float32 v_sequence = v_cvt_f32(vx_load(ptr_x_int)); + for (i = 0U; i < len_sub_vfloatStep; i+= v_float32_width) + { + v_float32 v_reg_i = v_add(vx_setall_f32((float)i), v_sequence); + + v_float32 v_x = v_sub(v_mul(v_x_start, v_reg_i, v_u), v_half); + + v_int32 v_x_int = v_trunc(v_x); + v_float32 v_x_fract = v_sub(v_x, v_cvt_f32(v_floor(v_x))); + + vx_store(ptr_x_int, v_x_int); + + for(j = 0U; j < v_float32_width; ++j) + ptr_for_cubicHermite[j] = *(tmp.ptr(0) + (ptr_x_int[j] - 1)); + v_float32 v_x_int_add_A = vx_load(ptr_for_cubicHermite); + + for(j = 0U; j < v_float32_width; ++j) + ptr_for_cubicHermite[j] = *(tmp.ptr(0) + (ptr_x_int[j])); + v_float32 v_x_int_add_B = vx_load(ptr_for_cubicHermite); + + for(j = 0U; j < v_float32_width; ++j) + ptr_for_cubicHermite[j] = *(tmp.ptr(0) + (ptr_x_int[j] + 1)); + v_float32 v_x_int_add_C = vx_load(ptr_for_cubicHermite); + + for(j = 0U; j < v_float32_width; ++j) + ptr_for_cubicHermite[j] = *(tmp.ptr(0) + (ptr_x_int[j] + 2)); + v_float32 v_x_int_add_D = vx_load(ptr_for_cubicHermite); + + + vx_store(dst.ptr(0) + i, simd_cubicHermite(v_x_int_add_A, v_x_int_add_B, v_x_int_add_C, v_x_int_add_D, v_x_fract)); + } +#endif + + // if the tail exists or scalar version + float *ptr = tmp.ptr(0) + 1U; + float lenScale = 1.0f / (float)(dstlen - 1U); + float U, X, xfract; + int xint; + for(; i < dstlen; ++i) + { + U = (float)i * lenScale; + X = (U * (float)srclen) - 0.5f; + xfract = X - floor(X); + xint = (int)X; + dst.at(i) = scal_cubicHermite(ptr[xint - 1], ptr[xint], ptr[xint + 1], ptr[xint + 2], xfract); + } + +} + +static void fir_f32(const float *pSrc, float *pDst, + const float *pCoeffs, float *pBuffer, + uint32_t numTaps, uint32_t blockSize) +{ + uint32_t copyLen = std::min(blockSize, numTaps); + + /////////////// delay line to the left /////////////// + uint32_t i = numTaps - 1U, k = 0U, j = 0U; + uint32_t value_i; + const float* ptr = pSrc + 1U - numTaps; + +#if (CV_SIMD || CV_SIMD_SCALABLE) + v_float32 v_pDst; + value_i = (uint32_t)std::max((int64_t)(numTaps + numTaps - 2U) - (int64_t)v_float32_width, (int64_t)0); + uint32_t value_k = (uint32_t)std::max((int64_t)copyLen - (int64_t)v_float32_width, (int64_t)0); + + + uint32_t value_j = (uint32_t)std::max((int64_t)(numTaps) - (int64_t)v_float32_width, (int64_t)0); + + for (; i < value_i && k < value_k; i += v_float32_width, k += v_float32_width) + { + v_float32 pSrc_data = vx_load(ptr + i); //vx_load(pSrc + (i + 1U - numTaps)); + vx_store(pBuffer + i, pSrc_data); + } +#endif + + // if the tail exists or scalar version + value_i = numTaps + numTaps - 2U; + for (; i < value_i && k < copyLen; ++i, ++k) + { + *(pBuffer + i) = *(ptr + i); // pBuffer[i] = pSrc[i + 1U - numTaps] + } + + + /////////////// process delay line /////////////// + i = 0U; k = 0U; + value_i = numTaps - 1U; + float *ptr_Buf; + + for(; i < value_i && k < copyLen; ++i, ++k) + { + ptr_Buf = pBuffer + i; + j = 0U; + +#if (CV_SIMD || CV_SIMD_SCALABLE) + + v_pDst = vx_setzero_f32(); + for (; j < value_j; j += v_float32_width) + { + v_float32 v_pBuffer = vx_load(ptr_Buf + j); //vx_load(pBuffer[i + j]) + v_float32 v_pCoeffs = vx_load(pCoeffs + j); //vx_load(pCoeffs[j]) + + v_pDst = v_fma(v_pBuffer, v_pCoeffs, v_pDst); // v_pDst = v_pBuffer * v_pCoeffs + v_pDst + } + pDst[i] = v_reduce_sum(v_pDst); +#endif + + // if the tail exists or scalar version + for (; j < numTaps; ++j) + pDst[i] += pCoeffs[j] * *(ptr_Buf + j); // pDst[i] += pCoeffs[j] * pBuffer[i + j]; + } + + + /////////////// process main block /////////////// + i = numTaps - 1U; + + for(; i < blockSize; ++i) + { + const float *ptr_Src = pSrc + (i + 1U - numTaps); + j = 0U; + +#if (CV_SIMD || CV_SIMD_SCALABLE) + v_pDst = vx_setzero_f32(); + for (; j < value_j; j += v_float32_width) + { + v_float32 v_pSrc = vx_load(ptr_Src + j); // vx_load(pSrc[i + j - (numTaps - 1)]) + v_float32 v_pCoeffs = vx_load(pCoeffs + j); //vx_load(pCoeffs[j]) + v_pDst = v_fma(v_pSrc, v_pCoeffs, v_pDst); + } + pDst[i] = v_reduce_sum(v_pDst); +#endif + + // if the tail exists or scalar version + for (; j < numTaps; ++j) + pDst[i] += pCoeffs[j] * *(ptr_Src + j); // pDst[i] += pCoeffs[j] * pSrc[i + j + 1U - numTaps]; + } + + + /////////////// move delay line left by copyLen elements /////////////// +#if (CV_SIMD || CV_SIMD_SCALABLE) + value_i = (uint32_t)std::max((int64_t)(numTaps - 1U) - (int64_t)v_float32_width, (int64_t)0); + ptr_Buf = pBuffer + copyLen; + + for(i = 0U; i < value_i; i += v_float32_width) + { + v_float32 v_pBuffer = vx_load(ptr_Buf + i); //vx_load(pBuffer[copyLen + i]) + vx_store(pBuffer + i, v_pBuffer); + } +#endif + + // if the tail exists or scalar version + value_i = numTaps - 1U; + for (; i < value_i; ++i) + { + pBuffer[i] = pBuffer[i + copyLen]; + } + + + /////////////// copy new elements /////////////// + /////////////// post-process delay line /////////////// + int l = (int)(numTaps - 2U); k = 0U; + +#if (CV_SIMD || CV_SIMD_SCALABLE) + int value_l = (int)v_float32_width; + const float* ptr_part = pSrc + (blockSize + 1U - numTaps - v_float32_width); + for(; l >= value_l && k < value_k; l -= value_l, k += v_float32_width) + { + v_float32 v_pSrc = vx_load(ptr_part + l); // vx_load(pSrc[blockSize - (numTaps - 1) + l - v_float32_width]) + vx_store(pBuffer + (l - value_l), v_pSrc); + } +#endif + const float* ptr_Src = pSrc + (blockSize + 1U - numTaps); + for(; l >= 0 && k < copyLen; --l, ++k) + { + pBuffer[l] = *(ptr_Src + l); // pBuffer[l] = pSrc[blockSize + 1U - numTaps + l]; + } +} + +void resampleSignal(InputArray inputSignal, OutputArray outputSignal, + const int inFreq, const int outFreq) +{ + CV_TRACE_FUNCTION(); + CV_Assert(!inputSignal.empty()); + CV_CheckGE(inFreq, 1000, ""); + CV_CheckGE(outFreq, 1000, ""); + if (inFreq == outFreq) + { + inputSignal.copyTo(outputSignal); + return; + } + uint32_t filtLen = 33U; + float beta = 3.395f; + std::vector filt_window(filtLen, 0.f); + init_filter(beta, filtLen, filt_window.data()); + float ratio = (float)outFreq / float(inFreq); + Mat1f inMat = inputSignal.getMat(); + Mat1f outMat = Mat1f(Size(cvFloor(inMat.cols * ratio), 1)); + cubicInterpolate(inMat, outMat.cols, outMat, inMat.cols); + if (inFreq < 2 * outFreq) + { + std::vector dlyl(filtLen * 2 - 1, 0.f); + std::vector ptmp(outMat.cols + 2 * filtLen, 0.); + + for (auto i = filtLen; i < outMat.cols + filtLen; ++i) + { + ptmp[i] = outMat.at(i - filtLen); + } + std::vector ptmp2(outMat.cols + 2 * filtLen, 0.f); + fir_f32(ptmp.data(), ptmp2.data(), filt_window.data(), dlyl.data(), filtLen, (uint32_t)(ptmp.size())); + for (auto i = filtLen; i < outMat.cols + filtLen; ++i) + { + outMat.at(i - filtLen) = ptmp2[i + cvFloor((float)filtLen / 2.f)]; + } + } + outputSignal.assign(std::move(outMat)); +} + + +} +} diff --git a/modules/signal/test/test_main.cpp b/modules/signal/test/test_main.cpp new file mode 100644 index 000000000..a6fc332d4 --- /dev/null +++ b/modules/signal/test/test_main.cpp @@ -0,0 +1,6 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. +#include "test_precomp.hpp" + +CV_TEST_MAIN("") diff --git a/modules/signal/test/test_precomp.hpp b/modules/signal/test/test_precomp.hpp new file mode 100644 index 000000000..c398e080f --- /dev/null +++ b/modules/signal/test/test_precomp.hpp @@ -0,0 +1,10 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. +#ifndef __OPENCV_TEST_PRECOMP_HPP__ +#define __OPENCV_TEST_PRECOMP_HPP__ + +#include "opencv2/ts.hpp" +#include "opencv2/signal.hpp" + +#endif diff --git a/modules/signal/test/test_signal.cpp b/modules/signal/test/test_signal.cpp new file mode 100644 index 000000000..377ac5734 --- /dev/null +++ b/modules/signal/test/test_signal.cpp @@ -0,0 +1,180 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. + +#include "test_precomp.hpp" +#include "opencv2/core/mat.hpp" +#include "opencv2/signal.hpp" + +#include +#include +#include +#include +#include +#include +#include + + +namespace opencv_test { namespace { + +using namespace cv; +using namespace cv::signal; + +float MSE(const Mat1f &outSignal, const Mat1f &refSignal) +{ + float mse = 0.f; + for (int i = 0; i < refSignal.cols; ++i) + { + mse += powf(outSignal.at(0,i) - refSignal.at(0,i), 2.f); + } + mse /= refSignal.cols; + return mse; +} + +// RRMSE = sqrt( MSE / SUM(sqr(refSignal(i))) ) * 100% +float RRMSE(const Mat1f &outSignal, const Mat1f &refSignal) +{ + float rrmse = 0.f; + float div = 0.f; + rrmse = MSE(outSignal, refSignal); + for (int i = 0; i < refSignal.cols; ++i) + { + div += powf(refSignal.at(0,i), 2.f); + } + rrmse /= div; + rrmse = sqrt(rrmse) * 100; + return rrmse; +} + +TEST(ResampleTest, simple_resample_test_up) +{ + Mat1f sample_signal(Size(1000U,1U)); + Mat1f outSignal; + std::iota(sample_signal.begin(), sample_signal.end(), 1.f); + resampleSignal(sample_signal, outSignal, 16000U, 32000U); + vector ref(outSignal.cols, 0.f); + for (uint32_t i = 0U; i < 2000U; ++i) + { + ref[i] = static_cast(i) / 2.f; + } + EXPECT_NEAR(cvtest::norm(ref, NORM_L2) / cvtest::norm(outSignal, NORM_L2), 1.0f, 0.05f) + << "\nL2_norm(refSignal) = " << cvtest::norm(ref, NORM_L2) + << "\nL2_norm(outSignal) = " << cvtest::norm(outSignal, NORM_L2); +} + +TEST(ResampleTest, resample_sin_signal_up_2) +{ + Mat1f sample_signal(Size(1000U,1U)); + Mat1f outSignal; + for (uint32_t i = 0U; i < (uint32_t)sample_signal.cols; ++i) + { + sample_signal.at(0, i) = sinf(float(i)); + } + resampleSignal(sample_signal, outSignal, 16000U, 32000U); + vector ref(outSignal.cols, 0.f); + for (uint32_t i = 0U; i < 2000U; ++i) + { + ref[i] = sin(static_cast(i) / 2.f); + } + EXPECT_NEAR(cvtest::norm(ref, NORM_L2) / cvtest::norm(outSignal, NORM_L2), 1.0f, 0.05f) + << "\nL2_norm(refSignal) = " << cvtest::norm(ref, NORM_L2) + << "\nL2_norm(outSignal) = " << cvtest::norm(outSignal, NORM_L2); +} + +TEST(ResampleTest, simple_resample_test_dn) +{ + Mat1f sample_signal(Size(1000U,1U)); + Mat1f outSignal; + std::iota(sample_signal.begin(), sample_signal.end(), 1.f); + resampleSignal(sample_signal, outSignal, 32000U, 16000U); + vector ref(outSignal.cols, 0.f); + for (uint32_t i = 0U; i < 500U; ++i) + { + ref[i] = static_cast(i) * 2.f; + } + EXPECT_NEAR(cvtest::norm(ref, NORM_L2) / cvtest::norm(outSignal, NORM_L2), 1.0f, 0.05f) + << "\nL2_norm(refSignal) = " << cvtest::norm(ref, NORM_L2) + << "\nL2_norm(outSignal) = " << cvtest::norm(outSignal, NORM_L2); +} + +TEST(ResampleTest, resample_sin_signal_dn_2) +{ + Mat1f sample_signal(Size(1000U,1U)); + Mat1f outSignal; + for (uint32_t i = 0U; i < (uint32_t)sample_signal.cols; ++i) + { + sample_signal.at(0, i) = sinf(float(i)); + } + resampleSignal(sample_signal, outSignal, 32000U, 16000U); + std::vector ref(outSignal.cols, 0.f); + for (uint32_t i = 0U; i < 500U; ++i) + { + ref[i] = sin(static_cast(i) * 2.f); + } + EXPECT_NEAR(cvtest::norm(ref, NORM_L2) / cvtest::norm(outSignal, NORM_L2), 1.0f, 0.05f) + << "\nL2_norm(refSignal) = " << cvtest::norm(ref, NORM_L2) + << "\nL2_norm(outSignal) = " << cvtest::norm(outSignal, NORM_L2); +} + +// produce 1s of signal @ freq hz +void fillSignal(uint32_t freq, Mat1f &inSignal) +{ + static std::default_random_engine e((unsigned int)(time(NULL))); + static std::uniform_real_distribution<> dis(0, 1); // range [0, 1) + static auto a = dis(e), b = dis(e), c = dis(e); + uint32_t idx = 0; + std::generate(inSignal.begin(), inSignal.end(), [&]() + { + float ret = static_cast(sin(idx/(float)freq + a) + 3 * sin(CV_PI / 4 * (idx/(float)freq + b)) + + 5 * sin(CV_PI/12 * idx/(float)freq + c) + 20*cos(idx/(float)freq*4000)); + idx++; + return ret; + }); +} + +class ResampleTestClass : public testing::TestWithParam> +{ +}; + +TEST_P(ResampleTestClass, func_test) { + auto params1 = GetParam(); + uint32_t inFreq = std::get<0>(params1); + uint32_t outFreq = std::get<1>(params1); + // 1 second @ inFreq hz + Mat1f inSignal(Size(inFreq, 1U)); + Mat1f outSignal; + // generating testing function as a sum of different sinusoids + fillSignal(inFreq, inSignal); + resampleSignal(inSignal, outSignal, inFreq, outFreq); + // reference signal + // 1 second @ outFreq hz + Mat1f refSignal(Size(outFreq, 1U)); + fillSignal(outFreq, refSignal); + // calculating maxDiff + float maxDiff = 0.f; + // exclude 2 elements and last 2 elements from testing + for (uint32_t i = 2; i < (uint32_t)refSignal.cols - 2; ++i) + { + if(maxDiff < abs(refSignal.at(0,i) - outSignal.at(0,i))) + { + maxDiff = abs(refSignal.at(0,i) - outSignal.at(0,i)); + } + } + auto max = std::max_element(outSignal.begin(), outSignal.end()); + float maxDiffRel = maxDiff / (*max); + EXPECT_LE(maxDiffRel, 0.35f); + // calculating relative error of L2 norms + EXPECT_NEAR(abs(cvtest::norm(outSignal, NORM_L2) - cvtest::norm(refSignal, NORM_L2)) / + cvtest::norm(refSignal, NORM_L2), 0.0f, 0.05f); + // calculating relative mean squared error + float rrmse = RRMSE(outSignal, refSignal); + // 1% error + EXPECT_LE(rrmse, 1.f); +} + +INSTANTIATE_TEST_CASE_P(RefSignalTestingCase, + ResampleTestClass, + ::testing::Combine(testing::Values(16000, 32000, 44100, 48000), + testing::Values(16000, 32000, 44100, 48000))); + +}} // namespace