amishutin
2 years ago
13 changed files with 696 additions and 0 deletions
@ -0,0 +1,2 @@ |
|||||||
|
set(the_description "Signal processing algorithms") |
||||||
|
ocv_define_module(signal opencv_core WRAP python) |
||||||
@ -0,0 +1,4 @@ |
|||||||
|
Signal Processing algorithms |
||||||
|
================================================ |
||||||
|
|
||||||
|
Signal resampling done with cubic interpolation and low-pass filtering |
||||||
@ -0,0 +1,17 @@ |
|||||||
|
// 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_HPP |
||||||
|
#define OPENCV_SIGNAL_HPP |
||||||
|
|
||||||
|
/**
|
||||||
|
* @defgroup signal Signal Processing |
||||||
|
* @{ |
||||||
|
* This module includes signal processing algorithms. |
||||||
|
* @} |
||||||
|
*/ |
||||||
|
|
||||||
|
#include "opencv2/core.hpp" |
||||||
|
#include "opencv2/signal/signal_resample.hpp" |
||||||
|
|
||||||
|
#endif |
||||||
@ -0,0 +1,32 @@ |
|||||||
|
// 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_SIGNAL_RESAMPLE_HPP |
||||||
|
#define OPENCV_SIGNAL_SIGNAL_RESAMPLE_HPP |
||||||
|
|
||||||
|
#include <opencv2/core.hpp> |
||||||
|
|
||||||
|
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 |
||||||
@ -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) |
||||||
@ -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 |
||||||
@ -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<uint32_t, uint32_t, uint32_t> > 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<float>(0, i) = sinf(float(i)); |
||||||
|
} |
||||||
|
declare.in(sample_signal).out(outSignal); |
||||||
|
TEST_CYCLE() resampleSignal(sample_signal, outSignal, inFreq, outFreq); |
||||||
|
SANITY_CHECK_NOTHING(); |
||||||
|
} |
||||||
|
|
||||||
|
}} |
||||||
@ -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 <opencv2/core.hpp> |
||||||
|
|
||||||
|
#endif |
||||||
@ -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 <opencv2/signal/signal_resample.hpp> |
||||||
|
#include <opencv2/core/mat.hpp> |
||||||
|
#include <opencv2/core/hal/intrin.hpp> |
||||||
|
#include <opencv2/core/utils/trace.hpp> |
||||||
|
|
||||||
|
#include <algorithm> |
||||||
|
#include <cmath> |
||||||
|
#include <vector> |
||||||
|
|
||||||
|
namespace cv { |
||||||
|
namespace signal { |
||||||
|
|
||||||
|
#if (CV_SIMD || CV_SIMD_SCALABLE) |
||||||
|
#define v_float32_width (uint32_t)VTraits<v_float32>::vlanes() |
||||||
|
const uint32_t v_float32_max_width = (uint32_t)VTraits<v_float32>::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<float>(i)); |
||||||
|
mul = mul / powf(tgammaf(static_cast<float>(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<float> 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<float>(0) = src.at<float>(0); |
||||||
|
|
||||||
|
#if (CV_SIMD || CV_SIMD_SCALABLE) |
||||||
|
v_float32 v_reg = vx_setall_f32(src.at<float>(srclen - 1U)); |
||||||
|
vx_store(tmp.ptr<float>(0) + (srclen - 1U), v_reg); |
||||||
|
#else // scalar version
|
||||||
|
tmp.at<float>(srclen + 1U) = src.at<float>(srclen - 1U); |
||||||
|
tmp.at<float>(srclen + 2U) = src.at<float>(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<float>(0) + i); |
||||||
|
vx_store(tmp.ptr<float>(0) + (i + 1U), v_copy); |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
// if the tail exists or scalar version
|
||||||
|
for (; i < srclen; ++i) |
||||||
|
{ |
||||||
|
tmp.at<float>(i + 1U) = src.at<float>(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<float>(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<float>(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<float>(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<float>(0) + (ptr_x_int[j] + 2)); |
||||||
|
v_float32 v_x_int_add_D = vx_load(ptr_for_cubicHermite); |
||||||
|
|
||||||
|
|
||||||
|
vx_store(dst.ptr<float>(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<float>(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<float>(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<float> 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<float> dlyl(filtLen * 2 - 1, 0.f); |
||||||
|
std::vector<float> ptmp(outMat.cols + 2 * filtLen, 0.); |
||||||
|
|
||||||
|
for (auto i = filtLen; i < outMat.cols + filtLen; ++i) |
||||||
|
{ |
||||||
|
ptmp[i] = outMat.at<float>(i - filtLen); |
||||||
|
} |
||||||
|
std::vector<float> 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<float>(i - filtLen) = ptmp2[i + cvFloor((float)filtLen / 2.f)]; |
||||||
|
} |
||||||
|
} |
||||||
|
outputSignal.assign(std::move(outMat)); |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
} |
||||||
|
} |
||||||
@ -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("") |
||||||
@ -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 |
||||||
@ -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 <vector> |
||||||
|
#include <numeric> |
||||||
|
#include <cmath> |
||||||
|
#include <string> |
||||||
|
#include <random> |
||||||
|
#include <ctime> |
||||||
|
#include <algorithm> |
||||||
|
|
||||||
|
|
||||||
|
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<float>(0,i) - refSignal.at<float>(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<float>(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<float> ref(outSignal.cols, 0.f); |
||||||
|
for (uint32_t i = 0U; i < 2000U; ++i) |
||||||
|
{ |
||||||
|
ref[i] = static_cast<float>(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<float>(0, i) = sinf(float(i)); |
||||||
|
} |
||||||
|
resampleSignal(sample_signal, outSignal, 16000U, 32000U); |
||||||
|
vector<float> ref(outSignal.cols, 0.f); |
||||||
|
for (uint32_t i = 0U; i < 2000U; ++i) |
||||||
|
{ |
||||||
|
ref[i] = sin(static_cast<float>(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<float> ref(outSignal.cols, 0.f); |
||||||
|
for (uint32_t i = 0U; i < 500U; ++i) |
||||||
|
{ |
||||||
|
ref[i] = static_cast<float>(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<float>(0, i) = sinf(float(i)); |
||||||
|
} |
||||||
|
resampleSignal(sample_signal, outSignal, 32000U, 16000U); |
||||||
|
std::vector<float> ref(outSignal.cols, 0.f); |
||||||
|
for (uint32_t i = 0U; i < 500U; ++i) |
||||||
|
{ |
||||||
|
ref[i] = sin(static_cast<float>(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<float>(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<std::tuple<int, int>> |
||||||
|
{ |
||||||
|
}; |
||||||
|
|
||||||
|
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<float>(0,i) - outSignal.at<float>(0,i))) |
||||||
|
{ |
||||||
|
maxDiff = abs(refSignal.at<float>(0,i) - outSignal.at<float>(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
|
||||||
Loading…
Reference in new issue