r8bbase.h

Go to the documentation of this file.

//$ nobt
 
#ifndef R8BBASE_INCLUDED
#define R8BBASE_INCLUDED
 
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include "r8bconf.h"
 
#if defined( _WIN32 )
    #include <windows.h>
#else // defined( _WIN32 )
    #include <pthread.h>
#endif // defined( _WIN32 )
 
#if defined( __SSE4_2__ ) || defined( __SSE4_1__ ) || \
    defined( __SSSE3__ ) || defined( __SSE3__ ) || defined( __SSE2__ ) || \
    defined( __x86_64__ ) || defined( __amd64 ) || defined( _M_X64 ) || \
    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
 
#elif defined( __aarch64__ ) || defined( __arm64 )
 
    #include <arm_neon.h>
 
    #define R8B_NEON
 
    #if !defined( __APPLE__ )
        #define R8B_SIMD_ISH // Shuffled interpolation is inefficient on M1.
    #endif // !defined( __APPLE__ )
 
#endif // AArch64
 
namespace r8b {
 
#define R8B_VERSION "6.5"
 
#define R8B_PI 3.14159265358979324
 
#define R8B_2PI 6.28318530717958648
 
#define R8B_3PI 9.42477796076937972
 
#define R8B_PId2 1.57079632679489662
 
#define R8BNOCTOR( ClassName ) \
    private: \
        ClassName( const ClassName& ) { } \
        ClassName& operator = ( const ClassName& ) { return( *this ); }
 
class CStdClassAllocator
{
public:
    void* operator new( const size_t, void* const p )
    {
        return( p );
    }
 
    void* operator new( const size_t n )
    {
        return( :: malloc( n ));
    }
 
    void* operator new[]( const size_t n )
    {
        return( :: malloc( n ));
    }
 
    void operator delete( void* const p )
    {
        :: free( p );
    }
 
    void operator delete[]( void* const p )
    {
        :: free( p );
    }
};
 
class CStdMemAllocator : public CStdClassAllocator
{
public:
    static void* allocmem( const size_t Size )
    {
        return( :: malloc( Size ));
    }
 
    static void* reallocmem( void* const p, const size_t Size )
    {
        return( :: realloc( p, Size ));
    }
 
    static void freemem( void* const p )
    {
        :: free( 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( NULL )
        , Data( NULL )
    {
    }
 
    CFixedBuffer( const int Capacity )
    {
        R8BASSERT( Capacity > 0 || Capacity == 0 );
 
        Data0 = allocmem( Capacity * sizeof( T ) + Alignment );
        Data = (T*) alignptr( Data0, Alignment );
 
        R8BASSERT( Data0 != NULL || Capacity == 0 );
    }
 
    ~CFixedBuffer()
    {
        freemem( Data0 );
    }
 
    void alloc( const int Capacity )
    {
        R8BASSERT( Capacity > 0 || Capacity == 0 );
 
        freemem( Data0 );
        Data0 = allocmem( Capacity * sizeof( T ) + Alignment );
        Data = (T*) alignptr( Data0, Alignment );
 
        R8BASSERT( Data0 != NULL || Capacity == 0 );
    }
 
    void realloc( const int PrevCapacity, const int NewCapacity )
    {
        R8BASSERT( PrevCapacity >= 0 );
        R8BASSERT( NewCapacity >= 0 );
 
        void* const NewData0 = allocmem( NewCapacity * sizeof( T ) +
            Alignment );
 
        T* const NewData = (T*) alignptr( NewData0, Alignment );
        const size_t CopySize = ( PrevCapacity > NewCapacity ?
            NewCapacity : PrevCapacity ) * sizeof( T );
 
        if( CopySize > 0 )
        {
            memcpy( NewData, Data, CopySize );
        }
 
        freemem( Data0 );
        Data0 = NewData0;
        Data = NewData;
 
        R8BASSERT( Data0 != NULL || NewCapacity == 0 );
    }
 
    void free()
    {
        freemem( Data0 );
        Data0 = NULL;
        Data = 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( 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 DelObj = Object;
        Object = NULL;
        delete DelObj;
    }
 
    T unkeep()
    {
        T ResObject = Object;
        Object = NULL;
        return( ResObject );
    }
 
private:
    T Object; 
};
 
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( NULL )
    {
    }
 
    CSyncKeeper( CSyncObject* const aSyncObj )
        : SyncObj( aSyncObj )
    {
        if( SyncObj != NULL )
        {
            SyncObj -> acquire();
        }
    }
 
    CSyncKeeper( CSyncObject& aSyncObj )
        : SyncObj( &aSyncObj )
    {
        SyncObj -> acquire();
    }
 
    ~CSyncKeeper()
    {
        if( SyncObj != NULL )
        {
            SyncObj -> release();
        }
    }
 
private:
    CSyncObject* SyncObj; 
};
 
#define R8BSYNC( SyncObject ) R8BSYNC_( SyncObject, __LINE__ )
#define R8BSYNC_( SyncObject, id ) R8BSYNC__( SyncObject, id )
#define R8BSYNC__( SyncObject, id ) 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 uint8_t 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 // min
 
#if !defined( max )
 
template< typename T >
inline T max( const T& v1, const T& v2 )
{
    return( v1 > v2 ? v1 : v2 );
}
 
#endif // max
 
inline double clampr( const double Value, const double minv,
    const double maxv )
{
    if( Value < minv )
    {
        return( minv );
    }
    else
    if( Value > maxv )
    {
        return( maxv );
    }
    else
    {
        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