r8bbase.h

Go to the documentation of this file.

//$ nobt

 
#ifndef R8BBASE_INCLUDED
#define R8BBASE_INCLUDED
 
#define R8B_VERSION "7.3" 


 
#if __cplusplus >= 201103L
 
    #include <cstdint>
    #include <cstring>
    #include <cmath>
    #include <mutex>
 
    #define R8B_CONST constexpr
    #define R8B_NULL nullptr
 
#else // __cplusplus >= 201103L
 
    #include <stdint.h>
    #include <string.h>
    #include <math.h>
 
    #if defined( _WIN32 )
        #include <Windows.h>
    #else // defined( _WIN32 )
        #include <pthread.h>
    #endif // defined( _WIN32 )
 
    #define R8B_CONST static const
    #define R8B_NULL NULL
 
#endif // __cplusplus >= 201103L
 
#include "r8bconf.h"
 
#if defined( __aarch64__ ) || defined( __arm64__ ) || \
    defined( _M_ARM64 ) || defined( _M_ARM64EC )
 
    #if defined( _MSC_VER )
        #include <arm64_neon.h>
    #else // defined( _MSC_VER )
        #include <arm_neon.h>
    #endif // defined( _MSC_VER )
 
    #define R8B_NEON
 
    #if !defined( __APPLE__ )
        #define R8B_SIMD_ISH // Shuffled interpolation is inefficient on M1.
    #endif // !defined( __APPLE__ )
 
#elif defined( __SSE2__ ) || defined( _M_AMD64 ) || \
    ( defined( _M_IX86_FP ) && _M_IX86_FP == 2 )
 
    #if defined( _MSC_VER )
        #include <intrin.h>
    #else // defined( _MSC_VER )
        #include <emmintrin.h>
    #endif // defined( _MSC_VER )
 
    #define R8B_SSE2
    #define R8B_SIMD_ISH
 
#endif // SSE2

 
#if defined( __clang__ )
 
    #define R8B_EXITDTOR __attribute__((always_destroy))
 
#else // defined( __clang__ )
 
    #define R8B_EXITDTOR
 
#endif // defined( __clang__ )

 
namespace r8b {
 
#if __cplusplus >= 201103L
 
    using std :: memcpy;
    using std :: memset;
    using std :: floor;
    using std :: ceil;
    using std :: exp;
    using std :: log;
    using std :: pow;
    using std :: fabs;
    using std :: sqrt;
    using std :: sin;
    using std :: cos;
    using std :: atan2;
    using std :: size_t;
    using std :: uintptr_t;
 
#endif // __cplusplus >= 201103L
 
R8B_CONST double R8B_PI = 3.14159265358979324; 
R8B_CONST double R8B_2PI = 6.28318530717958648; 
R8B_CONST double R8B_3PI = 9.42477796076937972; 
R8B_CONST double R8B_PId2 = 1.57079632679489662; 

 
#if __cplusplus >= 201103L
 
    #define R8BNOCTOR( ClassName ) \
        public: \
            ClassName( const ClassName& ) = delete; \
            ClassName& operator = ( const ClassName& ) = delete;
 
#else // __cplusplus >= 201103L
 
    #define R8BNOCTOR( ClassName ) \
        private: \
            ClassName( const ClassName& ) { } \
            ClassName& operator = ( const ClassName& ) { return( *this ); }
 
#endif // __cplusplus >= 201103L

 
class CStdClassAllocator
{
};

 
class CStdMemAllocator : public CStdClassAllocator
{
public:
 
    static void* allocmem( const size_t Size )
    {
        return( new char[ Size ]);
    }

 
    static void freemem( void* const p )
    {
        delete[] (char*) p;
    }
};

 
template< typename T >
inline T* alignptr( T* const ptr, const uintptr_t align )
{
    return( (T*) (( (uintptr_t) ptr + align - 1 ) & ~( align - 1 )));
}

 
template< typename T >
class CFixedBuffer : public R8B_MEMALLOCCLASS
{
    R8BNOCTOR( CFixedBuffer )
 
public:
    CFixedBuffer()
        : Data0( R8B_NULL )
        , Data( R8B_NULL )
    {
    }

 
    CFixedBuffer( const int Capacity )
    {
        R8BASSERT( Capacity >= 0 );
 
        Data0 = allocmem( (size_t) Capacity * sizeof( T ) + Alignment );
        Data = (T*) alignptr( Data0, Alignment );
 
        R8BASSERT( Data0 != R8B_NULL || Capacity == 0 );
    }
 
    ~CFixedBuffer()
    {
        freemem( Data0 );
    }

 
    void alloc( const int Capacity )
    {
        R8BASSERT( Capacity >= 0 );
 
        freemem( Data0 );
        Data0 = allocmem( (size_t) Capacity * sizeof( T ) + Alignment );
        Data = (T*) alignptr( Data0, Alignment );
 
        R8BASSERT( Data0 != R8B_NULL || Capacity == 0 );
    }

 
    void realloc( const int PrevCapacity, const int NewCapacity )
    {
        R8BASSERT( PrevCapacity >= 0 );
        R8BASSERT( NewCapacity >= 0 );
 
        void* const NewData0 = allocmem( (size_t) NewCapacity * sizeof( T ) +
            Alignment );
 
        T* const NewData = (T*) alignptr( NewData0, Alignment );
        const size_t CopySize = ( PrevCapacity > NewCapacity ?
            (size_t) NewCapacity : (size_t) PrevCapacity ) * sizeof( T );
 
        if( CopySize > 0 )
        {
            memcpy( NewData, Data, CopySize );
        }
 
        freemem( Data0 );
        Data0 = NewData0;
        Data = NewData;
 
        R8BASSERT( Data0 != R8B_NULL || NewCapacity == 0 );
    }

 
    void free()
    {
        freemem( Data0 );
        Data0 = R8B_NULL;
        Data = R8B_NULL;
    }

 
    operator T* () const
    {
        return( Data );
    }
 
private:
    static const size_t Alignment = 64; 
    void* Data0; 
    T* Data; 
};

 
template< typename T >
class CPtrKeeper
{
    R8BNOCTOR( CPtrKeeper )
 
public:
    CPtrKeeper()
        : Object( R8B_NULL )
    {
    }

 
    template< typename T2 >
    CPtrKeeper( T2 const aObject )
        : Object( aObject )
    {
    }
 
    ~CPtrKeeper()
    {
        delete Object;
    }

 
    template< typename T2 >
    void operator = ( T2 const aObject )
    {
        reset();
        Object = aObject;
    }

 
    T* operator -> () const
    {
        return( Object );
    }

 
    operator T* () const
    {
        return( Object );
    }

 
    void reset()
    {
        T* const DelObj = Object;
        Object = R8B_NULL;
        delete DelObj;
    }

 
    T* unkeep()
    {
        T* const ResObject = Object;
        Object = R8B_NULL;
        return( ResObject );
    }
 
private:
    T* Object; 
};
 
#if __cplusplus >= 201103L
 
typedef std :: mutex CSyncObject; 
typedef std :: lock_guard< std :: mutex > CSyncKeeper; 
 
#else // __cplusplus >= 201103L

 
class CSyncObject
{
    R8BNOCTOR( CSyncObject )
 
public:
    CSyncObject()
    {
        #if defined( _WIN32 )
            InitializeCriticalSectionAndSpinCount( &CritSec, 2000 );
        #else // defined( _WIN32 )
            pthread_mutexattr_t MutexAttrs;
            pthread_mutexattr_init( &MutexAttrs );
            pthread_mutexattr_settype( &MutexAttrs, PTHREAD_MUTEX_RECURSIVE );
            pthread_mutex_init( &Mutex, &MutexAttrs );
            pthread_mutexattr_destroy( &MutexAttrs );
        #endif // defined( _WIN32 )
    }
 
    ~CSyncObject()
    {
        #if defined( _WIN32 )
            DeleteCriticalSection( &CritSec );
        #else // defined( _WIN32 )
            pthread_mutex_destroy( &Mutex );
        #endif // defined( _WIN32 )
    }

 
    void acquire()
    {
        #if defined( _WIN32 )
            EnterCriticalSection( &CritSec );
        #else // defined( _WIN32 )
            pthread_mutex_lock( &Mutex );
        #endif // defined( _WIN32 )
    }

 
    void release()
    {
        #if defined( _WIN32 )
            LeaveCriticalSection( &CritSec );
        #else // defined( _WIN32 )
            pthread_mutex_unlock( &Mutex );
        #endif // defined( _WIN32 )
    }
 
private:
    #if defined( _WIN32 )
        CRITICAL_SECTION CritSec; 
    #else // defined( _WIN32 )
        pthread_mutex_t Mutex; 
    #endif // defined( _WIN32 )
};

 
class CSyncKeeper
{
    R8BNOCTOR( CSyncKeeper )
 
public:
    CSyncKeeper()
        : SyncObj( R8B_NULL )
    {
    }

 
    CSyncKeeper( CSyncObject* const aSyncObj )
        : SyncObj( aSyncObj )
    {
        if( SyncObj != R8B_NULL )
        {
            SyncObj -> acquire();
        }
    }

 
    CSyncKeeper( CSyncObject& aSyncObj )
        : SyncObj( &aSyncObj )
    {
        SyncObj -> acquire();
    }
 
    ~CSyncKeeper()
    {
        if( SyncObj != R8B_NULL )
        {
            SyncObj -> release();
        }
    }
 
private:
    CSyncObject* SyncObj; 
};
 
#endif // __cplusplus >= 201103L

 
#define R8BSYNC( SyncObject ) R8BSYNC1( SyncObject, __LINE__ )
#define R8BSYNC1( SyncObject, id ) R8BSYNC2( SyncObject, id ) 
#define R8BSYNC2( SyncObject, id ) \
    const CSyncKeeper SyncKeeper##id( SyncObject ) 

 
class CSineGen
{
public:
    CSineGen()
    {
    }

 
    CSineGen( const double si, const double ph )
        : svalue1( sin( ph ))
        , svalue2( sin( ph - si ))
        , sincr( 2.0 * cos( si ))
    {
    }

 
    CSineGen( const double si, const double ph, const double g )
        : svalue1( sin( ph ) * g )
        , svalue2( sin( ph - si ) * g )
        , sincr( 2.0 * cos( si ))
    {
    }

 
    void init( const double si, const double ph )
    {
        svalue1 = sin( ph );
        svalue2 = sin( ph - si );
        sincr = 2.0 * cos( si );
    }

 
    void init( const double si, const double ph, const double g )
    {
        svalue1 = sin( ph ) * g;
        svalue2 = sin( ph - si ) * g;
        sincr = 2.0 * cos( si );
    }

 
    double generate()
    {
        const double res = svalue1;
 
        svalue1 = sincr * res - svalue2;
        svalue2 = res;
 
        return( res );
    }
 
private:
    double svalue1; 
    double svalue2; 
    double sincr; 
};

 
inline int getBitOccupancy( const int v )
{
    static const unsigned char OccupancyTable[] =
    {
        1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
        5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
        6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
        6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
        7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
        7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
        7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
        7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
        8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
        8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
        8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
        8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
        8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
        8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
        8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
        8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8
    };
 
    const int tt = v >> 16;
 
    if( tt != 0 )
    {
        const int t = v >> 24;
 
        return( t != 0 ? 24 + OccupancyTable[ t & 0xFF ] :
            16 + OccupancyTable[ tt ]);
    }
    else
    {
        const int t = v >> 8;
 
        return( t != 0 ? 8 + OccupancyTable[ t ] : OccupancyTable[ v ]);
    }
}

 
inline void calcFIRFilterResponse( const double* flt, int fltlen,
    const double th, double& re0, double& im0, const int fltlat = 0 )
{
    const double sincr = 2.0 * cos( th );
    double cvalue1;
    double svalue1;
 
    if( fltlat == 0 )
    {
        cvalue1 = 1.0;
        svalue1 = 0.0;
    }
    else
    {
        cvalue1 = cos( -fltlat * th );
        svalue1 = sin( -fltlat * th );
    }
 
    double cvalue2 = cos( -( fltlat + 1 ) * th );
    double svalue2 = sin( -( fltlat + 1 ) * th );
 
    double re = 0.0;
    double im = 0.0;
 
    while( fltlen > 0 )
    {
        re += cvalue1 * flt[ 0 ];
        im += svalue1 * flt[ 0 ];
        flt++;
        fltlen--;
 
        double tmp = cvalue1;
        cvalue1 = sincr * cvalue1 - cvalue2;
        cvalue2 = tmp;
 
        tmp = svalue1;
        svalue1 = sincr * svalue1 - svalue2;
        svalue2 = tmp;
    }
 
    re0 = re;
    im0 = im;
}

 
inline double calcFIRFilterGroupDelay( const double* const flt,
    const int fltlen, const double th )
{
    const int Count = 2;
    const double thd2 = 1e-9;
    double ths[ Count ] = { th - thd2, th + thd2 }; // Side-band frequencies.
 
    if( ths[ 0 ] < 0.0 )
    {
        ths[ 0 ] = 0.0;
    }
 
    if( ths[ 1 ] > R8B_PI )
    {
        ths[ 1 ] = R8B_PI;
    }
 
    double ph1[ Count ];
    int i;
 
    for( i = 0; i < Count; i++ )
    {
        double re1;
        double im1;
 
        calcFIRFilterResponse( flt, fltlen, ths[ i ], re1, im1 );
        ph1[ i ] = atan2( im1, re1 );
    }
 
    if( fabs( ph1[ 1 ] - ph1[ 0 ]) > R8B_PI )
    {
        if( ph1[ 1 ] > ph1[ 0 ])
        {
            ph1[ 1 ] -= R8B_2PI;
        }
        else
        {
            ph1[ 1 ] += R8B_2PI;
        }
    }
 
    const double thd = ths[ 1 ] - ths[ 0 ];
 
    return(( ph1[ 1 ] - ph1[ 0 ]) / thd );
}

 
inline void normalizeFIRFilter( double* const p, const int l,
    const double DCGain, const int pstep = 1 )
{
    R8BASSERT( l > 0 );
    R8BASSERT( pstep != 0 );
 
    double s = 0.0;
    double* pp = p;
    int i = l;
 
    while( i > 0 )
    {
        s += *pp;
        pp += pstep;
        i--;
    }
 
    s = DCGain / s;
    pp = p;
    i = l;
 
    while( i > 0 )
    {
        *pp *= s;
        pp += pstep;
        i--;
    }
}

 
inline void calcSpline3p8Coeffs( double* const c, const double xm3,
    const double xm2, const double xm1, const double x0, const double x1,
    const double x2, const double x3, const double x4 )
{
    c[ 0 ] = x0;
    c[ 1 ] = ( 61.0 * ( x1 - xm1 ) + 16.0 * ( xm2 - x2 ) +
        3.0 * ( x3 - xm3 )) * 1.31578947368421052e-2;
 
    c[ 2 ] = ( 106.0 * ( xm1 + x1 ) + 10.0 * x3 + 6.0 * xm3 - 3.0 * x4 -
        29.0 * ( xm2 + x2 ) - 167.0 * x0 ) * 1.31578947368421052e-2;
 
    c[ 3 ] = ( 91.0 * ( x0 - x1 ) + 45.0 * ( x2 - xm1 ) +
        13.0 * ( xm2 - x3 ) + 3.0 * ( x4 - xm3 )) * 1.31578947368421052e-2;
}

 
inline void calcSpline2p8Coeffs( double* const c, const double xm3,
    const double xm2, const double xm1, const double x0, const double x1,
    const double x2, const double x3, const double x4 )
{
    c[ 0 ] = x0;
    c[ 1 ] = ( 61.0 * ( x1 - xm1 ) + 16.0 * ( xm2 - x2 ) +
        3.0 * ( x3 - xm3 )) * 1.31578947368421052e-2;
 
    c[ 2 ] = ( 106.0 * ( xm1 + x1 ) + 10.0 * x3 + 6.0 * xm3 - 3.0 * x4 -
        29.0 * ( xm2 + x2 ) - 167.0 * x0 ) * 1.31578947368421052e-2;
}

 
inline void calcSpline3p4Coeffs( double* const c, const double* const y )
{
    c[ 0 ] = y[ 1 ];
    c[ 1 ] = 0.5 * ( y[ 2 ] - y[ 0 ]);
    c[ 2 ] = y[ 0 ] - 2.5 * y[ 1 ] + y[ 2 ] + y[ 2 ] - 0.5 * y[ 3 ];
    c[ 3 ] = 0.5 * ( y[ 3 ] - y[ 0 ] ) + 1.5 * ( y[ 1 ] - y[ 2 ]);
}

 
inline void calcSpline3p6Coeffs( double* const c, const double* const y )
{
    c[ 0 ] = y[ 2 ];
    c[ 1 ] = ( 11.0 * ( y[ 3 ] - y[ 1 ]) + 2.0 * ( y[ 0 ] - y[ 4 ])) / 14.0;
    c[ 2 ] = ( 20.0 * ( y[ 1 ] + y[ 3 ]) + 2.0 * y[ 5 ] - 4.0 * y[ 0 ] -
        7.0 * y[ 4 ] - 31.0 * y[ 2 ]) / 14.0;
 
    c[ 3 ] = ( 17.0 * ( y[ 2 ] - y[ 3 ]) + 9.0 * ( y[ 4 ] - y[ 1 ]) +
        2.0 * ( y[ 0 ] - y[ 5 ])) / 14.0;
}
 
#if !defined( min )

 
template< typename T >
inline T min( const T& v1, const T& v2 )
{
    return( v1 < v2 ? v1 : v2 );
}
 
#endif // !defined( min )
 
#if !defined( max )

 
template< typename T >
inline T max( const T& v1, const T& v2 )
{
    return( v1 > v2 ? v1 : v2 );
}
 
#endif // !defined( max )

 
inline double clampr( const double Value, const double minv,
    const double maxv )
{
    if( Value < minv )
    {
        return( minv );
    }
 
    if( Value > maxv )
    {
        return( maxv );
    }
 
    return( Value );
}

 
inline double sqr( const double x )
{
    return( x * x );
}

 
inline double pow_a( const double v, const double p )
{
    return( exp( p * log( fabs( v ) + 1e-300 )));
}

 
inline double gauss( const double v )
{
    return( exp( -( v * v )));
}

 
inline double asinh( const double v )
{
    return( log( v + sqrt( v * v + 1.0 )));
}

 
inline double besselI0( const double x )
{
    const double ax = fabs( x );
    double y;
 
    if( ax < 3.75 )
    {
        y = x / 3.75;
        y *= y;
 
        return( 1.0 + y * ( 3.5156229 + y * ( 3.0899424 + y * ( 1.2067492 +
            y * ( 0.2659732 + y * ( 0.360768e-1 + y * 0.45813e-2 ))))));
    }
 
    y = 3.75 / ax;
 
    return( exp( ax ) / sqrt( ax ) * ( 0.39894228 + y * ( 0.1328592e-1 +
        y * ( 0.225319e-2 + y * ( -0.157565e-2 + y * ( 0.916281e-2 +
        y * ( -0.2057706e-1 + y * ( 0.2635537e-1 + y * ( -0.1647633e-1 +
        y * 0.392377e-2 )))))))));
}
 
} // namespace r8b
 
#endif // R8BBASE_INCLUDED