From 9546f1dd83409c45f55915527dad5e59a046eebc Mon Sep 17 00:00:00 2001 From: camilo Date: Mon, 15 Jan 2024 12:33:23 -0500 Subject: [PATCH] added interp1 class. bump to 1.0.8 --- check/main.cpp | 19 ++ check/qlibs_cpp_test.cpp | 107 ++++++++- check/sa_check.ewp | 9 + check/sa_check.ewt | 9 + keywords.txt | 13 +- library.json | 2 +- library.properties | 2 +- src/CMakeLists.txt | 1 + src/fis.cpp | 3 +- src/include/interp1.hpp | 152 +++++++++++++ src/interp1.cpp | 470 +++++++++++++++++++++++++++++++++++++++ src/qlibs.h | 7 +- src/smoother.cpp | 24 +- 13 files changed, 802 insertions(+), 16 deletions(-) create mode 100644 src/include/interp1.hpp create mode 100644 src/interp1.cpp diff --git a/check/main.cpp b/check/main.cpp index 1ff684c..68e345e 100644 --- a/check/main.cpp +++ b/check/main.cpp @@ -3,5 +3,24 @@ using namespace qlibs; int main() { + + mat<2,2> x( + 1.0f, 2.0f, + 3.0f, 4.0f + ); + mat<2,1> y( + 1.0f, + 3.0f + ); + mat<2,1> z( + 2.0f, + 2.0f + ); + + mat<2,1> result = 4.5 - z - 2.0f*x*y - 3; + x*=x; + result.display(); + x.display(); + return 0; } diff --git a/check/qlibs_cpp_test.cpp b/check/qlibs_cpp_test.cpp index 25c7915..e2d8ed8 100644 --- a/check/qlibs_cpp_test.cpp +++ b/check/qlibs_cpp_test.cpp @@ -11,7 +11,8 @@ void test_fp16( void ); void test_crc( void ); void test_ltisys( void ); void test_ffmath( void ); - +void test_mat( void ); +void test_interp1( void ); void test_fis3( void ) { @@ -304,6 +305,108 @@ void test_ffmath(void) cout << ffmath::getInf() << endl; } +void test_mat( void ) +{ + cout<<"MAT TEST"< x( + 1.0f, 2.0f, + 3.0f, 4.0f + ); + mat<2,1> y( + 1.0f, + 3.0f + ); + mat<2,1> z( + 2.0f, + 2.0f + ); + + mat<2,1> result = 4.5 - z - 2.0f*x*y - 3; + mat<4,4> I( MAT_IDENTITY ); + x*=x; + result.display(); + x.display(); + cout<< x(2) << endl; + auto j = !y; + j.display(); + I.display(); + auto ix = x.inv(); + x.display(); + ix.display(); + (x*ix).display(); + +} + +void test_interp1( void ) +{ + real_t tx[] = { 1.0f, 2.0f, 3.0f, 4.0f }; + real_t ty[] = { 5.0f, 9.0f, 12.0f, 15.0f }; + interp1 interpolation( tx, ty ); + + cout << "linear" << endl; + interpolation.setMethod( INTERP1_LINEAR ); + cout << interpolation.get( 2.5 ) << endl; + cout << interpolation.get( 3.1 ) << endl; + cout << interpolation.get( 0.5 ) << endl; + cout << interpolation.get( 5.0 ) << endl; + + cout << "sine" << endl; + interpolation.setMethod( INTERP1_SINE ); + cout << interpolation.get( 2.5 ) << endl; + cout << interpolation.get( 3.1 ) << endl; + cout << interpolation.get( 0.5 ) << endl; + cout << interpolation.get( 5.0 ) << endl; + + cout << "cubic" << endl; + interpolation.setMethod( INTERP1_CUBIC ); + cout << interpolation.get( 2.5 ) << endl; + cout << interpolation.get( 3.1 ) << endl; + cout << interpolation.get( 0.5 ) << endl; + cout << interpolation.get( 5.0 ) << endl; + + cout << "hermite" << endl; + interpolation.setMethod( INTERP1_HERMITE ); + cout << interpolation.get( 2.5 ) << endl; + cout << interpolation.get( 3.1 ) << endl; + cout << interpolation.get( 0.5 ) << endl; + cout << interpolation.get( 5.0 ) << endl; + + cout << "nearest" << endl; + interpolation.setMethod( INTERP1_NEAREST); + cout << interpolation.get( 2.5 ) << endl; + cout << interpolation.get( 3.1 ) << endl; + cout << interpolation.get( 0.5 ) << endl; + cout << interpolation.get( 5.0 ) << endl; + + cout << "next" << endl; + interpolation.setMethod( INTERP1_NEXT); + cout << interpolation.get( 2.5 ) << endl; + cout << interpolation.get( 3.1 ) << endl; + cout << interpolation.get( 0.5 ) << endl; + cout << interpolation.get( 5.0 ) << endl; + + cout << "previous" << endl; + interpolation.setMethod( INTERP1_PREVIOUS); + cout << interpolation.get( 2.5 ) << endl; + cout << interpolation.get( 3.1 ) << endl; + cout << interpolation.get( 0.5 ) << endl; + cout << interpolation.get( 5.0 ) << endl; + + cout << "spline" << endl; + interpolation.setMethod( INTERP1_SPLINE); + cout << interpolation.get( 2.5 ) << endl; + cout << interpolation.get( 3.1 ) << endl; + cout << interpolation.get( 0.5 ) << endl; + cout << interpolation.get( 5.0 ) << endl; + + cout << "cspline" << endl; + interpolation.setMethod( INTERP1_CONSTRAINED_SPLINE); + cout << interpolation.get( 2.5 ) << endl; + cout << interpolation.get( 3.1 ) << endl; + cout << interpolation.get( 0.5 ) << endl; + cout << interpolation.get( 5.0 ) << endl; +} + int main() { test_crc(); @@ -314,5 +417,7 @@ int main() test_fis2(); test_fis3(); test_ffmath(); + test_mat(); + test_interp1(); return 0; } \ No newline at end of file diff --git a/check/sa_check.ewp b/check/sa_check.ewp index 81f6c3c..c36d1e2 100644 --- a/check/sa_check.ewp +++ b/check/sa_check.ewp @@ -2164,9 +2164,15 @@ $PROJ_DIR$\..\src\include\generic.hpp + + $PROJ_DIR$\..\src\include\interp1.hpp + $PROJ_DIR$\..\src\include\ltisys.hpp + + $PROJ_DIR$\..\src\include\mat.hpp + $PROJ_DIR$\..\src\include\mathex.hpp @@ -2210,6 +2216,9 @@ $PROJ_DIR$\..\src\generic.cpp + + $PROJ_DIR$\..\src\interp1.cpp + $PROJ_DIR$\..\src\ltisys.cpp diff --git a/check/sa_check.ewt b/check/sa_check.ewt index 76427a3..9d0becb 100644 --- a/check/sa_check.ewt +++ b/check/sa_check.ewt @@ -2868,9 +2868,15 @@ $PROJ_DIR$\..\src\include\generic.hpp + + $PROJ_DIR$\..\src\include\interp1.hpp + $PROJ_DIR$\..\src\include\ltisys.hpp + + $PROJ_DIR$\..\src\include\mat.hpp + $PROJ_DIR$\..\src\include\mathex.hpp @@ -2914,6 +2920,9 @@ $PROJ_DIR$\..\src\generic.cpp + + $PROJ_DIR$\..\src\interp1.cpp + $PROJ_DIR$\..\src\ltisys.cpp diff --git a/keywords.txt b/keywords.txt index ba51a0d..a0376ab 100644 --- a/keywords.txt +++ b/keywords.txt @@ -6,6 +6,8 @@ # Datatypes(KEYWORD1) ####################################### +interp1 KEYWORD1 +interp1Method KEYWORD1 real_t KEYWORD1 bitfield KEYWORD1 crc KEYWORD1 @@ -65,6 +67,7 @@ tdl KEYWORD1 # Methods and Functions(KEYWORD2) ####################################### +setMethod KEYWORD2 bitfield_size KEYWORD2 setup KEYWORD2 clearAll KEYWORD2 @@ -224,7 +227,15 @@ insertSample KEYWORD2 ####################################### # Constants(LITERAL1) - +INTERP1_NEXT LITERAL1 +INTERP1_PREVIOUS LITERAL1 +INTERP1_NEAREST LITERAL1 +INTERP1_LINEAR LITERAL1 +INTERP1_SINE LITERAL1 +INTERP1_CUBIC LITERAL1 +INTERP1_HERMITE LITERAL1 +INTERP1_SPLINE LITERAL1 +INTERP1_CONSTRAINED_SPLINE LITERAL1 CRC8 LITERAL1 CRC16 LITERAL1 CRC32 LITERAL1 diff --git a/library.json b/library.json index 6753f1b..ef7c568 100644 --- a/library.json +++ b/library.json @@ -16,7 +16,7 @@ "maintainer": true } ], - "version": "1.0.7", + "version": "1.0.8", "license": "MIT", "frameworks": "arduino", "platforms": "*" diff --git a/library.properties b/library.properties index b55547d..a718405 100644 --- a/library.properties +++ b/library.properties @@ -1,5 +1,5 @@ name=qlibs -version=1.0.7 +version=1.0.8 license=MIT author=J. Camilo Gomez C. maintainer=J. Camilo Gomez C. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b474fcf..ec74ca3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,6 +21,7 @@ target_sources( ${PROJECT_NAME} INTERFACE "${CMAKE_CURRENT_LIST_DIR}/bitfield.cpp" "${CMAKE_CURRENT_LIST_DIR}/generic.cpp" "${CMAKE_CURRENT_LIST_DIR}/ffmath.cpp" + "${CMAKE_CURRENT_LIST_DIR}/interp1.cpp" ) target_compile_options( ${PROJECT_NAME} INTERFACE diff --git a/src/fis.cpp b/src/fis.cpp index 1110d10..1567c61 100644 --- a/src/fis.cpp +++ b/src/fis.cpp @@ -145,8 +145,9 @@ bool fis::instance::setupOutput( const fis::tag t, /*cstat -CERT-STR34-C*/ xOutput[ t ].min = Min; xOutput[ t ].max = Max; + /*cstat -CERT-FLP36-C*/ xOutput[ t ].res = ( xOutput[ t ].max - xOutput[ t ].min )/static_cast( nPoints ); - /*cstat +CERT-STR34-C*/ + /*cstat +CERT-STR34-C +CERT-FLP36-C*/ retVal = true; } diff --git a/src/include/interp1.hpp b/src/include/interp1.hpp new file mode 100644 index 0000000..a4ffc38 --- /dev/null +++ b/src/include/interp1.hpp @@ -0,0 +1,152 @@ +/*! + * @file interp1.hpp + * @author J. Camilo Gomez C. + * @version 1.01 + * @note This file is part of the qLibs-cpp distribution. + * @brief Helper class that implements a Tapped Delay Line (TDL). A TDL is a + * delay line that provides access to its contents at arbitrary intermediate + * delay length values. + * This class runs in constant time O(1), so it becomes useful when you need to + * work with long delayed lines. + **/ + +#ifndef QLIBS_INTERP1 +#define QLIBS_INTERP1 + +#include + +/** +* @brief The qLibs++ library namespace. +*/ +namespace qlibs { + /** @addtogroup qinterp1 1D Interpolation + * @brief One-dimensional interpolation class. + * @{ + */ + + /** + * @brief An enum with all the available interpolation methods. + */ + enum interp1Method { + INTERP1_NEXT = 0, /*!< Return the next neighbor.*/ + INTERP1_PREVIOUS, /*!< Return the previous neighbor.*/ + INTERP1_NEAREST, /*!< Return the nearest neighbor.*/ + INTERP1_LINEAR, /*!< Linear interpolation from nearest neighbors.*/ + INTERP1_SINE, /*!< Sine interpolation.*/ + INTERP1_CUBIC, /*!< Cubic interpolation.*/ + INTERP1_HERMITE, /*!< Piecewise cubic Hermite interpolation.*/ + INTERP1_SPLINE, /*!< Catmull spline interpolation.*/ + INTERP1_CONSTRAINED_SPLINE, /*!< A special kind of spline that doesn't overshoot.*/ + INTERP1_MAX, + }; + + /** + * @brief A 1D interpolation object. + */ + class interp1 { + private: + using interp1Fcn_t = real_t (*)( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + const real_t *xData{ nullptr }; + const real_t *yData{ nullptr }; + const size_t dataSize{ 0U }; + interp1Fcn_t method{ &linear }; + static real_t slope( const real_t * const tx, + const real_t * const ty, + const size_t i ); + static real_t firstDerivate( const real_t * const tx, + const real_t * const ty, + const size_t n, + const size_t i ); + static real_t leftSecondDerivate( const real_t * const tx, + const real_t * const ty, + const size_t n, + const size_t i ); + static real_t rightSecondDerivate( const real_t * const tx, + const real_t * const ty, + const size_t n, + const size_t i ); + static real_t next( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + static real_t previous( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + static real_t nearest( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + static real_t linear( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + static real_t sine( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + static real_t cubic( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + static real_t hermite( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + static real_t spline( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + static real_t cSpline( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ); + public: + virtual ~interp1() {} + /** + * @brief Constructor for the 1D interpolation instance. + * @param[in] xTable An array of size @a sizeTable with the x points sorted in ascending order. + * @param[in] yTabLe An array of size @a sizeTable with the y points. + * @param[in] sizeTable The number of points in @a xTable @a yTable + */ + interp1( const real_t * const xTable, + const real_t * const yTabLe, + const size_t sizeTable ) : xData( xTable ), yData( yTabLe ), dataSize( sizeTable ) {} + + /** + * @brief Constructor for the 1D interpolation instance. + * @param[in] xTable An array of size @a sizeTable with the x points sorted in ascending order. + * @param[in] yTabLe An array of size @a sizeTable with the y points. + */ + template + interp1( real_t (&xTable)[ sizeTable ], + real_t (&yTable)[ sizeTable ] ) : interp1( xTable, yTable, sizeTable ) {} + + /** + * @brief Specify the interpolation method to use. + * @param[in] m The interpolation method. + * @return @c true on success otherwise @c false. + */ + bool setMethod( const interp1Method m ); + template + + /** + * @brief Interpolate input point @a x to determine the value of @a y + * at the points @a xi using the current method. If value if beyond + * the endpoints, extrapolation is performed using the current method. + * @return @c The interpolated-extrapolated @a y value. + */ + inline real_t get( const T x ) + { + return method( static_cast( x ), xData, yData, dataSize ); + } + }; + + /** @}*/ +} + + +#endif /*QLIBS_INTERP1*/ \ No newline at end of file diff --git a/src/interp1.cpp b/src/interp1.cpp new file mode 100644 index 0000000..09123e8 --- /dev/null +++ b/src/interp1.cpp @@ -0,0 +1,470 @@ +#include +#include +#include + +using namespace qlibs; + +/*cstat -CERT-INT30-C_a*/ +/*============================================================================*/ +bool interp1::setMethod( const interp1Method m ) +{ + bool retValue = false; + static const interp1Fcn_t im[ ] = { + &interp1::next, + &interp1::previous, + &interp1::nearest, + &interp1::linear, + &interp1::sine, + &interp1::cubic, + &interp1::hermite, + &interp1::spline, + &interp1::cSpline, + }; + + if ( ( m > 0 ) && ( m < INTERP1_MAX ) ) { + method = im[ m ]; + retValue = true; + } + + return retValue; +} +/*============================================================================*/ +real_t interp1::next( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ) +{ + real_t y = ffmath::getNan(); + + if ( ( tableSize >= 2U ) && ( nullptr != tx ) && ( nullptr != ty ) ) { + size_t nearestIndex = tableSize - 1U; + + for ( size_t i = 0U ; i < tableSize - 1U ; ++i ) { + if ( x < tx[ i + 1U ] ) { + nearestIndex = i; + break; + } + } + + if ( x >= tx[ tableSize - 1U ] ) { + y = ty[ tableSize - 1U ]; + } + else { + y = ty[ nearestIndex + 1U ]; + } + } + + return y; +} +/*============================================================================*/ +real_t interp1::previous( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ) +{ + real_t y = ffmath::getNan(); + + if ( ( tableSize >= 2U ) && ( nullptr != tx ) && ( nullptr != ty ) ) { + size_t nearestIndex = 0U; + + if ( x <= tx[ 0 ] ) { + y = ty[ 0 ]; + } + else { + for ( size_t i = 1U ; i < tableSize; ++i ) { + if ( x < tx [ i ] ) { + break; + } + nearestIndex = i; + } + y = ty[ nearestIndex ]; + } + } + return y; +} +/*============================================================================*/ +real_t interp1::nearest( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ) +{ + real_t y = ffmath::getNan(); + + if ( ( tableSize >= 2U ) && ( nullptr != tx ) && ( nullptr != ty ) ) { + size_t nearestIndex = 0U; + real_t minDistance = ffmath::absf( x - tx[ 0 ] ); + + for ( size_t i = 1U ; i < tableSize; ++i ) { + const real_t distance = ffmath::absf( x - tx[ i ] ); + + if ( distance <= minDistance) { + minDistance = distance; + nearestIndex = i; + } + } + y = ty[ nearestIndex ]; + } + return y; +} +/*============================================================================*/ +real_t interp1::linear( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ) +{ + real_t y = ffmath::getNan(); + + if ( ( tableSize >= 2U ) && ( nullptr != tx ) && ( nullptr != ty ) ) { + if ( x < tx[ 0 ] ) { + const real_t x0 = tx[ 0 ]; + const real_t x1 = tx[ 1 ]; + const real_t y0 = ty[ 0 ]; + const real_t y1 = ty[ 1 ]; + y = y0 + ( ( y1 - y0 )/( x1 - x0 ) )*( x - x0 ); + } + else if ( x > tx[ tableSize - 1U ] ) { + const real_t x0 = tx[ tableSize - 2U ]; + const real_t x1 = tx[ tableSize - 1U ]; + const real_t y0 = ty[ tableSize - 2U ]; + const real_t y1 = ty[ tableSize - 1U ]; + y = y1 + ( ( y0 - y1 )/( x0 - x1 ) )*( x - x1 ); + } + else { + const int maxIndex = static_cast( tableSize ) - 1; + for ( int i = 0; i < maxIndex; ++i ) { + if ( ( x >= tx[ i ] ) && ( x <= tx[ i + 1 ] ) ) { + const real_t x0 = tx[ i ]; + const real_t x1 = tx[ i + 1 ]; + const real_t y0 = ty[ i ]; + const real_t y1 = ty[ i + 1 ]; + y = y0 + ( ( y1 - y0 )/( x1 - x0 ) )*( x - x0 ); + break; + } + } + } + } + return y; +} +/*============================================================================*/ +real_t interp1::sine( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ) +{ + real_t y = ffmath::getNan(); + if ( ( tableSize >= 2U ) && ( nullptr != tx ) && ( nullptr != ty ) ) { + if ( x < tx[ 0 ] ) { + const real_t x0 = tx[ 0 ]; + const real_t x1 = tx[ 1 ]; + const real_t y0 = ty[ 0 ]; + const real_t y1 = ty[ 1 ]; + const real_t w = 0.5_re - 0.5_re*ffmath::cos( ffmath::FFP_PI*( x - x0 )/( x1 - x0 ) ); + y = y0 + w*( y1 - y0 ); + } + else if ( x > tx[ tableSize - 1U ] ) { + const real_t x0 = tx[ tableSize - 2U ]; + const real_t x1 = tx[ tableSize - 1U ]; + const real_t y0 = ty[ tableSize - 2U ]; + const real_t y1 = ty[ tableSize - 1U ]; + const real_t w = 0.5_re - 0.5_re*ffmath::cos( ffmath::FFP_PI*( x - x1 )/( x0 - x1 ) ); + y = y1 + w*( y0 - y1 ); + } + else { + for ( size_t i = 1; i < tableSize; ++i ) { + if ( x <= tx[ i ] ) { + const real_t x0 = tx[ i - 1 ]; + const real_t x1 = tx[ i ]; + const real_t y0 = ty[ i - 1 ]; + const real_t y1 = ty[ i]; + const real_t w = 0.5_re - 0.5_re*ffmath::cos( ffmath::FFP_PI*( x - x0 )/( x1 - x0 ) ); + y = y0 + w*( y1 - y0 ); + } + } + } + } + return y; +} +/*============================================================================*/ +real_t interp1::cubic( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ) +{ + real_t y = ffmath::getNan(); + if ( ( tableSize >= 4U ) && ( nullptr != tx ) && ( nullptr != ty ) ) { + if ( x < tx[ 0 ] ) { + const real_t x0 = tx[ 0 ]; + const real_t x1 = tx[ 1 ]; + const real_t y0 = ty[ 0 ]; + const real_t y1 = ty[ 1 ]; + const real_t h = x1 - x0; + const real_t t = ( x - x0 )/h; + const real_t t2 = t*t; + const real_t t3 = t2*t; + + y = ( 2.0_re*t3 - 3.0_re*t2 + 1.0_re )*y0 + + ( t3 - 2.0_re*t2 + t )*h*( y0 - y1 ) + + ( -2.0_re*t3 + 3.0_re*t2 )*y1 + + ( t3 - t2 )*h*( y1 - y0 ); + } + else if ( x > tx[ tableSize - 1U ] ) { + const real_t x0 = tx[ tableSize - 2U ]; + const real_t x1 = tx[ tableSize - 1U ]; + const real_t y0 = ty[ tableSize - 2U ]; + const real_t y1 = ty[ tableSize - 1U ]; + const real_t h = x1 - x0; + const real_t t = ( x - x1 )/h; + const real_t t2 = t*t; + const real_t t3 = t2*t; + + y = ( 2.0_re*t3 - 3.0_re*t2 + 1.0_re )*y1 + + ( t3 - 2.0_re*t2 + t )*h*( y0 - ty[ tableSize - 3U ] ) + + ( -2.0_re*t3 + 3.0_re*t2 )*y0 + + ( t3 - t2 )*h*( y1 - y0 ); + } + else { + for ( size_t i = 1U; i < tableSize; ++i ) { + if ( x <= tx[ i ] ) { + const real_t x0 = tx[ i - 1U ]; + const real_t x1 = tx[ i ]; + const real_t y0 = ty[ i - 1U ]; + const real_t y1 = ty[ i ]; + const real_t h = x1 - x0; + const real_t t = ( x - x0 )/h; + const real_t t2 = t*t; + const real_t t3 = t2*t; + y = ( 2.0_re*t3 - 3.0_re*t2 + 1.0_re )*y0 + + ( t3 - 2.0_re*t2 + t )*h*( y0 - ty[ i - 2U ] ) + + ( -2.0_re*t3 + 3.0_re*t2 )*y1 + + ( t3 - t2 )*h*( y1 - y0 ); + break; + } + } + } + } + return y; +} +/*============================================================================*/ +real_t interp1::hermite( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ) +{ + real_t y = ffmath::getNan(); + + if ( ( tableSize >= 2U ) && ( nullptr != tx ) && ( nullptr != ty ) ) { + if ( x < tx[ 0 ] ) { + const real_t x0 = tx[ 0 ]; + const real_t x1 = tx[ 1 ]; + const real_t y0 = ty[ 0 ]; + const real_t y1 = ty[ 1 ]; + y = y0 + ( ( y1 - y0 )/( x1 - x0 ) )*( x - x0 ); + } + else if ( x > tx[ tableSize -1U ] ) { + const real_t x0 = tx[ tableSize - 2U ]; + const real_t x1 = tx[ tableSize - 1U ]; + const real_t y0 = ty[ tableSize - 2U ]; + const real_t y1 = ty[ tableSize - 1U ]; + y = y1 + ( ( y0 - y1 )/( x0 - x1 ) )*( x - x1 ); + } + else { + y = 0.0_re; + for ( size_t i = 0U ; i < tableSize; ++i ) { + real_t term = ty[ i ]; + + for ( size_t j = 0U ; j < tableSize; ++j ) { + if ( i != j ) { + term *= ( x - tx[ j ] )/( tx[ i ] - tx[ j ] ); + } + } + + y += term; + } + } + } + return y; +} +/*============================================================================*/ +real_t interp1::spline( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ) +{ + real_t y = ffmath::getNan(); + + if ( ( tableSize >= 4U ) && ( nullptr != tx ) && ( nullptr != ty ) ) { + size_t i = 0U; + /* Extrapolation for x beyond the range*/ + if ( x <= tx[ 0 ] ) { + i = 0U; + } + else if ( x >= tx[ tableSize - 1U ] ) { + i = tableSize - 2U; /* Use the last interval for extrapolation*/ + } + else { + while ( x >= tx[ i + 1U ] ) { + i++; + } + } + + if ( isEqual( x , tx[ i + 1U ] ) ) { + y = ty[ i + 1U ]; + } + else { + const real_t t = ( x - tx[ i ] )/( tx[ i + 1U ] - tx[ i ] ); + const real_t t_2 = t*t; + const real_t t_3 = t_2*t; + const real_t h00 = 2.0_re*t_3 - 3.0_re*t_2 + 1.0_re; + const real_t h10 = t_3 - 2.0_re*t_2 + t; + const real_t h01 = 3.0_re * t_2 - 2.0_re*t_3; + const real_t h11 = t_3 - t_2; + const real_t x1_x0 = tx[ i + 1U ] - tx[ i ]; + const real_t y0 = ty[ i ]; + const real_t y1 = ty[ i + 1 ]; + real_t m0, m1; + + if ( 0U == i ) { + m0 = ( ty[ 1 ] - ty[ 0 ] )/( tx[ 1 ] - tx[ 0 ] ); + m1 = ( ty[ 2 ] - ty[ 0 ] )/( tx[ 2 ] - tx[ 0 ] ); + } + else if ( ( tableSize - 2U ) == i ) { + m0 = ( ty[ tableSize - 1U ] - ty[ tableSize - 3U ] )/( tx[ tableSize - 1U ] - tx[ tableSize - 3U ] ); + m1 = ( ty[ tableSize - 1U ] - ty[ tableSize - 2U ] )/( tx[ tableSize - 1U ] - tx[ tableSize - 2U ] ); + } + else { + m0 = slope( tx, ty, i ); + m1 = slope( tx, ty, i + 1U ); + } + y = ( h00*y0 ) + ( h01*y1 ) + ( h10*x1_x0*m0 ) + ( h11*x1_x0*m1 ); + } + } + + return y; +} +/*============================================================================*/ +real_t interp1::cSpline( const real_t x, + const real_t * const tx, + const real_t * const ty, + const size_t tableSize ) +{ + real_t y = ffmath::getNan(); + + if ( ( tableSize >= 4U ) && ( nullptr != tx ) && ( nullptr != ty ) ) { + size_t i = 0U; + /* Extrapolation for x beyond the range*/ + if ( x <= tx[ 0 ] ) { + i = 0U; + } + else if ( x >= tx[ tableSize - 1U ] ) { + i = tableSize - 2U; /* Use the last interval for extrapolation*/ + } + else { + while ( x >= tx[ i + 1U ] ) { + i++; + } + } + + if ( isEqual( x , tx[ i + 1U ] ) ) { + y = ty[ i + 1U ]; + } + else { + const real_t x0 = tx[ i ]; + const real_t x1 = tx[ i + 1U ]; + const real_t y0 = ty[ i ]; + const real_t y1 = ty[ i + 1U ]; + const real_t fd2i_xl1 = leftSecondDerivate( tx, ty, tableSize - 1, i + 1 ); + const real_t fd2i_x = rightSecondDerivate( tx, ty, tableSize - 1, i + 1 ); + const real_t x0_x1 = x0 - x1; + const real_t x1_2 = x1*x1; + const real_t x1_3 = x1_2*x1; + const real_t x0_2 = x0*x0; + const real_t x0_3 = x0_2*x0; + const real_t d = ( fd2i_x - fd2i_xl1 )/( 6.0_re*x0_x1 ); + const real_t c = ( ( x0*fd2i_xl1 ) - ( x1*fd2i_x ) )/( 2.0_re*x0_x1 ); + const real_t b = ( y0 - y1 - c*( x0_2 - x1_2 ) - d*( x0_3 - x1_3 ) )/x0_x1; + const real_t a = y1 - ( b*x1 ) - ( c*x1_2 ) - ( d*x1_3 ); + y = a + x*( b + x*( c + ( x*d ) ) ); + } + } + + return y; +} +/*============================================================================*/ +real_t interp1::slope( const real_t * const tx, + const real_t * const ty, + const size_t i ) +{ + real_t m; + + if ( isEqual( tx[ i + 1U ], tx[ i - 1U ] ) ) { + m = 0.0_re; + } + else { + m = ( ty[ i + 1U ] - ty[ i - 1U ] )/( tx[ i + 1U ] - tx[ i - 1U ] ); + } + + return m; +} +/*============================================================================*/ +real_t interp1::firstDerivate( const real_t * const tx, + const real_t * const ty, + const size_t n, + const size_t i ) +{ + real_t fd1_x; + + if ( 0 == i ) { + const real_t dx = tx[ 1 ] - tx[ 0 ]; + const real_t dy = ty[ 1 ] - ty[ 0 ]; + fd1_x = 1.5_re*( dy/dx ); + fd1_x -= 1.0_re/( ( tx[ 2 ] - tx[ 0 ] )/( ty[ 2 ] - ty[ 0 ] ) + dx/dy ); + } + else if ( n == i ) { + const real_t dx = tx[ n ] - tx[ n - 1 ]; + const real_t dy = ty[ n ] - ty[ n - 1 ]; + fd1_x = 1.5_re*( dy/dx ); + fd1_x -= 1.0_re/( ( tx[ n ] - tx[ n - 2 ] )/( ty[ n ] - ty[ n - 2 ] ) + ( dx/dy ) ); + } + else { + const real_t tmp1 = ( tx[ i + 1 ] - tx[ i ] )/( ty[ i + 1 ] - ty[ i ] ); + const real_t tmp2 = ( tx[ i ] - tx[ i - 1 ] )/( ty[ i ] - ty[ i - 1 ] ); + + if ( ( tmp1*tmp2 ) < 0.0_re ) { + fd1_x = 0.0_re; + } + else{ + fd1_x = 2.0_re/( tmp1 + tmp2 ); + } + } + return fd1_x; +} +/*============================================================================*/ +real_t interp1::leftSecondDerivate( const real_t * const tx, + const real_t * const ty, + const size_t n, + const size_t i ) +{ + const real_t fdi_x = firstDerivate( tx, ty, n, i ); + const real_t fdi_xl1 = firstDerivate( tx, ty, n, i - 1U ); + const real_t xi_delta = tx[ i ] - tx[ i - 1U ]; + real_t fd2l_x = -2.0_re*( fdi_x + ( 2.0_re*fdi_xl1 ) )/xi_delta; + + fd2l_x += 6.0_re*( ty[ i ] - ty[ i - 1U ] )/( xi_delta*xi_delta ); + return fd2l_x; +} +/*============================================================================*/ +real_t interp1::rightSecondDerivate( const real_t * const tx, + const real_t * const ty, + const size_t n, + const size_t i ) +{ + const real_t fdi_x = firstDerivate( tx, ty, n, i ); + const real_t fdi_xl1 = firstDerivate( tx, ty, n, i - 1U ); + const real_t xi_delta = tx[ i ] - tx[ i - 1U ]; + real_t fd2r_x = 2.0_re*( ( 2.0_re*fdi_x ) + fdi_xl1 )/xi_delta; + + fd2r_x -= 6.0_re*( ty[ i ] - ty[ i - 1U ] )/( xi_delta*xi_delta ); + return fd2r_x; +} +/*============================================================================*/ +/*cstat +CERT-INT30-C_a*/ diff --git a/src/qlibs.h b/src/qlibs.h index e20903a..64df885 100644 --- a/src/qlibs.h +++ b/src/qlibs.h @@ -41,8 +41,8 @@ This file is part of the QuarkTS++ OS distribution. #ifndef QLIBS_CPP_H #define QLIBS_CPP_H -#define QLIBS_CPP_VERSION "1.0.7" -#define QLIBS_CPP_VERNUM ( 107U ) +#define QLIBS_CPP_VERSION "1.0.8" +#define QLIBS_CPP_VERNUM ( 108U ) #define QLIBS_CPP_CAPTION "qLibs++" QLIBS_CPP_VERSION #include @@ -59,7 +59,8 @@ This file is part of the QuarkTS++ OS distribution. #include #include #include - +#include +#include using namespace qlibs; diff --git a/src/smoother.cpp b/src/smoother.cpp index 34ce014..5205c55 100644 --- a/src/smoother.cpp +++ b/src/smoother.cpp @@ -102,8 +102,9 @@ real_t smootherMWM1::smooth( const real_t x ) windowSet( w, wsize, x ); init = false; } - + /*cstat -CERT-FLP36-C*/ return discreteSystem::updateFIR( w, wsize, x )/static_cast( wsize ); + /*cstat +CERT-FLP36-C*/ } /*============================================================================*/ bool smootherMWM2::setup( real_t *window, @@ -121,8 +122,9 @@ bool smootherMWM2::setup( real_t *window, /*============================================================================*/ real_t smootherMWM2::smooth( const real_t x ) { + /*cstat -CERT-FLP36-C*/ const real_t wsize = static_cast( itemCount ); - + /*cstat +CERT-FLP36-C*/ if ( init ) { flush( x ); sum = x*wsize; @@ -165,7 +167,9 @@ real_t smootherMOR1::smooth( const real_t x ) w[ 0 ] = m; /*replace the outlier with the dynamic median*/ } /*compute new mean for next iteration*/ + /*cstat -CERT-FLP36-C*/ m = ( mc + w[ 0 ] ) / static_cast( wsize ); + /*cstat +CERT-FLP36-C*/ return w[ 0 ]; } /*============================================================================*/ @@ -189,8 +193,9 @@ bool smootherMOR2::setup( real_t *window, real_t smootherMOR2::smooth( const real_t x ) { real_t xx = x; + /*cstat -CERT-FLP36-C*/ const real_t wsize = static_cast( itemCount ); - + /*cstat +CERT-FLP36-C*/ if ( init ) { flush( x ); sum = wsize*x; @@ -216,20 +221,20 @@ bool smootherGMWF::setup( const real_t sg, { bool retValue = false; const size_t ws = wk_size; - + /*cstat -CERT-FLP36-C*/ if ( ( nullptr != window ) && ( wk_size > 0U ) && ( c < static_cast( ws ) ) && ( sg > 0.0_re ) ) { real_t r, sum = 0.0_re; size_t i; real_t l, center; /*cstat -MISRAC++2008-5-0-7*/ l = static_cast( wk_size - 1U ); - /*cstat +MISRAC++2008-5-0-7*/ + /*cstat +MISRAC++2008-5-0-7 +CERT-FLP36-C*/ center = c - l; r = 2.0_re*sg*sg; for ( i = 0U ; i < ws ; ++i ) { - /*cstat -MISRAC++2008-5-0-7*/ + /*cstat -MISRAC++2008-5-0-7 -CERT-FLP36-C*/ real_t d = static_cast( i ) - l; /*symmetry*/ - /*cstat +MISRAC++2008-5-0-7*/ + /*cstat +MISRAC++2008-5-0-7 +CERT-FLP36-C*/ d -= center; /*cstat -CERT-FLP32-C_b*/ kernel[ i ] = ffmath::exp( -( d*d )/r ); @@ -334,7 +339,9 @@ bool smootherDESF::setup( const real_t a, if ( ( a > 0.0_re ) && ( a < 1.0_re ) && ( b > 0.0_re ) && ( b < 1.0_re ) ) { alpha = a; beta = b; + /*cstat -CERT-FLP36-C */ n = static_cast( nS ); + /*cstat +CERT-FLP36-C */ retValue = reset(); } @@ -382,8 +389,9 @@ real_t smootherALNF::smooth( const real_t x ) real_t xe; if ( init ) { + /*cstat -CERT-FLP36-C */ const real_t np = 1.0_re/static_cast( n ); - + /*cstat +CERT-FLP36-C */ windowSet( xx, n, x ); windowSet( w, n, np ); if ( nullptr != w_1 ) {