CDSPSincFilterGen.h

Go to the documentation of this file.

//$ nobt
//$ nocpp

 
#ifndef R8B_CDSPSINCFILTERGEN_INCLUDED
#define R8B_CDSPSINCFILTERGEN_INCLUDED
 
#include "r8bbase.h"
 
namespace r8b {

 
class CDSPSincFilterGen
{
public:
    double Len2; 
    double Len2i; 
    int KernelLen; 
    int fl2; 
 
    double Freq1; 
    double Freq2; 
    double FracDelay; 

 
    enum EWindowFunctionType
    {
        wftCosine, 
        wftKaiser, 
        wftGaussian 
    };
 
    typedef double( CDSPSincFilterGen ::* CWindowFunc )(); 

 
    void initWindow( const EWindowFunctionType WinType = wftCosine,
        const double* const Params = R8B_NULL, const bool UsePower = false )
    {
        R8BASSERT( Len2 >= 2.0 );
 
        fl2 = (int) floor( Len2 );
        KernelLen = fl2 + fl2 + 1;
 
        setWindow( WinType, Params, UsePower, true );
    }

 
    void initBand( const EWindowFunctionType WinType = wftCosine,
        const double* const Params = R8B_NULL, const bool UsePower = false )
    {
        R8BASSERT( Len2 >= 2.0 );
 
        fl2 = (int) floor( Len2 );
        KernelLen = fl2 + fl2 + 1;
 
        setWindow( WinType, Params, UsePower, true );
    }

 
    void initHilbert( const EWindowFunctionType WinType = wftCosine,
        const double* const Params = R8B_NULL, const bool UsePower = false )
    {
        R8BASSERT( Len2 >= 2.0 );
 
        fl2 = (int) floor( Len2 );
        KernelLen = fl2 + fl2 + 1;
 
        setWindow( WinType, Params, UsePower, true );
    }

 
    void initFrac( const EWindowFunctionType WinType = wftCosine,
        const double* const Params = R8B_NULL, const bool UsePower = false )
    {
        R8BASSERT( Len2 >= 2.0 );
 
        fl2 = (int) ceil( Len2 );
        KernelLen = fl2 + fl2;
 
        setWindow( WinType, Params, UsePower, false, FracDelay );
    }

 
    double calcWindowHann()
    {
        return( 0.5 + 0.5 * w1.generate() );
    }

 
    double calcWindowHamming()
    {
        return( 0.54 + 0.46 * w1.generate() );
    }

 
    double calcWindowBlackman()
    {
        return( 0.42 + 0.5 * w1.generate() + 0.08 * w2.generate() );
    }

 
    double calcWindowNuttall()
    {
        return( 0.355768 + 0.487396 * w1.generate() +
            0.144232 * w2.generate() + 0.012604 * w3.generate() );
    }

 
    double calcWindowBlackmanNuttall()
    {
        return( 0.3635819 + 0.4891775 * w1.generate() +
            0.1365995 * w2.generate() + 0.0106411 * w3.generate() );
    }

 
    double calcWindowKaiser()
    {
        const double n = 1.0 - sqr( wn * Len2i + KaiserLen2Frac );
        wn++;
 
        if( n <= 0.0 )
        {
            return( 0.0 );
        }
 
        return( besselI0( KaiserBeta * sqrt( n )) * KaiserMul );
    }

 
    double calcWindowGaussian()
    {
        const double f = exp( -0.5 * sqr( wn * GaussianSigmaI +
            GaussianSigmaFrac ));
 
        wn++;
 
        return( f );
    }

 
    void generateWindow( double* op,
        CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman )
    {
        op += fl2;
        double* op2 = op;
 
        int l = fl2;
 
        if( Power < 0.0 )
        {
            *op = ( *this.*wfunc )();
 
            while( l > 0 )
            {
                const double v = ( *this.*wfunc )();
 
                op++;
                op2--;
                *op = v;
                *op2 = v;
                l--;
            }
        }
        else
        {
            *op = pow_a(( *this.*wfunc )(), Power );
 
            while( l > 0 )
            {
                const double v = pow_a(( *this.*wfunc )(), Power );
 
                op++;
                op2--;
                *op = v;
                *op2 = v;
                l--;
            }
        }
    }

 
    void generateBand( double* op,
        CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman )
    {
        CSineGen f2( Freq2, 0.0, 1.0 / R8B_PI );
        f2.generate();
 
        op += fl2;
        double* op2 = op;
        const double pw = Power;
        int t = 1;
 
        if( Freq1 < 2.3e-13 )
        {
            if( pw < 0.0 )
            {
                *op = Freq2 * ( *this.*wfunc )() / R8B_PI;
 
                while( t <= fl2 )
                {
                    const double v = f2.generate() * ( *this.*wfunc )() / t;
                    op++;
                    op2--;
                    *op = v;
                    *op2 = v;
                    t++;
                }
            }
            else
            {
                *op = Freq2 * pow_a(( *this.*wfunc )(), pw ) / R8B_PI;
 
                while( t <= fl2 )
                {
                    const double v = f2.generate() *
                        pow_a(( *this.*wfunc )(), pw ) / t;
 
                    op++;
                    op2--;
                    *op = v;
                    *op2 = v;
                    t++;
                }
            }
        }
        else
        {
            CSineGen f1( Freq1, 0.0, 1.0 / R8B_PI );
            f1.generate();
 
            if( pw < 0.0 )
            {
                *op = ( Freq2 - Freq1 ) * ( *this.*wfunc )() / R8B_PI;
 
                while( t <= fl2 )
                {
                    const double v = ( f2.generate() - f1.generate() ) *
                        ( *this.*wfunc )() / t;
 
                    op++;
                    op2--;
                    *op = v;
                    *op2 = v;
                    t++;
                }
            }
            else
            {
                *op = ( Freq2 - Freq1 ) *
                    pow_a(( *this.*wfunc )(), pw ) / R8B_PI;
 
                while( t <= fl2 )
                {
                    const double v = ( f2.generate() - f1.generate() ) *
                        pow_a(( *this.*wfunc )(), pw ) / t;
 
                    op++;
                    op2--;
                    *op = v;
                    *op2 = v;
                    t++;
                }
            }
        }
    }

 
    void generateHilbert( double* op,
        CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman )
    {
        static const double fvalues[ 2 ] = { 0.0, 2.0 / R8B_PI };
        op += fl2;
        double* op2 = op;
 
        ( *this.*wfunc )();
        *op = 0.0;
 
        int t = 1;
 
        if( Power < 0.0 )
        {
            while( t <= fl2 )
            {
                const double v = fvalues[ t & 1 ] * ( *this.*wfunc )() / t;
                op++;
                op2--;
                *op = v;
                *op2 = -v;
                t++;
            }
        }
        else
        {
            while( t <= fl2 )
            {
                const double v = fvalues[ t & 1 ] *
                    pow_a(( *this.*wfunc )(), Power ) / t;
 
                op++;
                op2--;
                *op = v;
                *op2 = -v;
                t++;
            }
        }
    }

 
    void generateFrac( double* op,
        CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman,
        const int opinc = 1 )
    {
        R8BASSERT( opinc != 0 );
 
        const double pw = Power;
        const double fd = FracDelay;
        int t = -fl2;
 
        if( t + fd < -Len2 )
        {
            ( *this.*wfunc )();
            *op = 0.0;
            op += opinc;
            t++;
        }
 
        double f = sin( fd * R8B_PI ) / R8B_PI;
 
        if(( t & 1 ) != 0 )
        {
            f = -f;
        }
 
        int IsZeroX = ( fabs( fd - 1.0 ) < 2.3e-13 );
        int mt = 0 - IsZeroX;
        IsZeroX = ( IsZeroX || fabs( fd ) < 2.3e-13 );
 
        if( pw < 0.0 )
        {
            while( t < mt )
            {
                *op = f * ( *this.*wfunc )() / ( t + fd );
                op += opinc;
                t++;
                f = -f;
            }
 
            if( IsZeroX ) // t+FracDelay==0
            {
                *op = ( *this.*wfunc )();
            }
            else
            {
                *op = f * ( *this.*wfunc )() / fd; // t==0
            }
 
            mt = fl2 - 2;
 
            while( t < mt )
            {
                op += opinc;
                t++;
                f = -f;
                *op = f * ( *this.*wfunc )() / ( t + fd );
            }
 
            op += opinc;
            t++;
            f = -f;
            const double ut = t + fd;
            *op = ( ut > Len2 ? 0.0 : f * ( *this.*wfunc )() / ut );
        }
        else
        {
            while( t < mt )
            {
                *op = f * pow_a(( *this.*wfunc )(), pw ) / ( t + fd );
                op += opinc;
                t++;
                f = -f;
            }
 
            if( IsZeroX ) // t+FracDelay==0
            {
                *op = pow_a(( *this.*wfunc )(), pw );
            }
            else
            {
                *op = f * pow_a(( *this.*wfunc )(), pw ) / fd; // t==0
            }
 
            mt = fl2 - 2;
 
            while( t < mt )
            {
                op += opinc;
                t++;
                f = -f;
                *op = f * pow_a(( *this.*wfunc )(), pw ) / ( t + fd );
            }
 
            op += opinc;
            t++;
            f = -f;
            const double ut = t + FracDelay;
            *op = ( ut > Len2 ? 0.0 :
                f * pow_a(( *this.*wfunc )(), pw ) / ut );
        }
    }
 
private:
    double Power; 
    CSineGen w1; 
    CSineGen w2; 
    CSineGen w3; 
 
    double KaiserBeta; 
    double KaiserMul; 
    double KaiserLen2Frac; 
    double GaussianSigmaI; 
    double GaussianSigmaFrac; 
 
    int wn; 

 
    void setWindowKaiser( const double* Params, const bool UsePower,
        const bool IsCentered )
    {
        wn = ( IsCentered ? 0 : -fl2 );
 
        if( Params == R8B_NULL )
        {
            KaiserBeta = 9.5945013206755156;
            Power = ( UsePower ? 1.9718457932433306 : -1.0 );
        }
        else
        {
            KaiserBeta = clampr( Params[ 0 ], 1.0, 350.0 );
            Power = ( UsePower ? fabs( Params[ 1 ]) : -1.0 );
        }
 
        KaiserMul = 1.0 / besselI0( KaiserBeta );
        Len2i = 1.0 / Len2;
        KaiserLen2Frac = FracDelay * Len2i;
    }

 
    void setWindowGaussian( const double* Params, const bool UsePower,
        const bool IsCentered )
    {
        wn = ( IsCentered ? 0 : -fl2 );
 
        if( Params == R8B_NULL )
        {
            GaussianSigmaI = 1.0;
            Power = -1.0;
        }
        else
        {
            GaussianSigmaI = clampr( fabs( Params[ 0 ]), 1e-1, 100.0 );
            Power = ( UsePower ? fabs( Params[ 1 ]) : -1.0 );
        }
 
        GaussianSigmaI *= Len2;
        GaussianSigmaI = 1.0 / GaussianSigmaI;
        GaussianSigmaFrac = FracDelay * GaussianSigmaI;
    }

 
    void setWindow( const EWindowFunctionType WinType,
        const double* const Params, const bool UsePower,
        const bool IsCentered, const double UseFracDelay = 0.0 )
    {
        FracDelay = UseFracDelay;
 
        if( WinType == wftCosine )
        {
            if( IsCentered )
            {
                w1.init( R8B_PI / Len2, R8B_PId2 );
                w2.init( R8B_2PI / Len2, R8B_PId2 );
                w3.init( R8B_3PI / Len2, R8B_PId2 );
            }
            else
            {
                const double step1 = R8B_PI / Len2;
                w1.init( step1, R8B_PId2 - step1 * fl2 + step1 * FracDelay );
 
                const double step2 = R8B_2PI / Len2;
                w2.init( step2, R8B_PId2 - step2 * fl2 + step2 * FracDelay );
 
                const double step3 = R8B_3PI / Len2;
                w3.init( step3, R8B_PId2 - step3 * fl2 + step3 * FracDelay );
            }
 
            Power = ( UsePower && Params != R8B_NULL ? Params[ 0 ] : -1.0 );
        }
        else
        if( WinType == wftKaiser )
        {
            setWindowKaiser( Params, UsePower, IsCentered );
        }
        else
        if( WinType == wftGaussian )
        {
            setWindowGaussian( Params, UsePower, IsCentered );
        }
    }
};
 
} // namespace r8b
 
#endif // R8B_CDSPSINCFILTERGEN_INCLUDED