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; 
    int KernelLen; 
    int fl2; 

    union
    {
        struct
        {
            double Freq1; 
            double Freq2; 
        };

        struct
        {
            double FracDelay; 
        };
    };

    enum EWindowFunctionType
    {
        wftCosine, 
        wftKaiser, 
        wftGaussian 
    };

    typedef double( CDSPSincFilterGen :: *CWindowFunc )(); 

    void initWindow( const EWindowFunctionType WinType = wftCosine,
        const double* const Params = 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 = NULL, const bool UsePower = false )
    {
        R8BASSERT( Len2 >= 2.0 );

        fl2 = (int) floor( Len2 );
        KernelLen = fl2 + fl2 + 1;
        f1.init( Freq1, 0.0 );
        f2.init( Freq2, 0.0 );

        setWindow( WinType, Params, UsePower, true );
    }

    void initHilbert( const EWindowFunctionType WinType = wftCosine,
        const double* const Params = 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 = 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 / Len2 + KaiserLen2Frac );
        wn++;

        if( n < 0.0 )
        {
            return( 0.0 );
        }

        return( besselI0( KaiserBeta * sqrt( n )) / KaiserDiv );
    }

    double calcWindowGaussian()
    {
        const double f = exp( -0.5 * sqr( wn / GaussianSigma +
            GaussianSigmaFrac ));

        wn++;

        return( f );
    }

    template< class T >
    void generateWindow( T* op,
        CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman )
    {
        op += fl2;
        T* 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 = pows(( *this.*wfunc )(), Power );

            while( l > 0 )
            {
                const double v = pows(( *this.*wfunc )(), Power );

                op++;
                op2--;
                *op = v;
                *op2 = v;
                l--;
            }
        }
    }

    template< class T >
    void generateBand( T* op,
        CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman )
    {
        op += fl2;
        T* op2 = op;
        f1.generate();
        f2.generate();
        int t = 1;

        if( Power < 0.0 )
        {
            *op = ( Freq2 - Freq1 ) * ( *this.*wfunc )() / M_PI;

            while( t <= fl2 )
            {
                const double v = ( f2.generate() - f1.generate() ) *
                    ( *this.*wfunc )() / t / M_PI;

                op++;
                op2--;
                *op = v;
                *op2 = v;
                t++;
            }
        }
        else
        {
            *op = ( Freq2 - Freq1 ) * pows(( *this.*wfunc )(), Power ) / M_PI;

            while( t <= fl2 )
            {
                const double v = ( f2.generate() - f1.generate() ) *
                    pows(( *this.*wfunc )(), Power ) / t / M_PI;

                op++;
                op2--;
                *op = v;
                *op2 = v;
                t++;
            }
        }
    }

    template< class T >
    void generateHilbert( T* op,
        CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman )
    {
        static const double fvalues[ 2 ] = { 0.0, 2.0 };
        op += fl2;
        T* 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 / M_PI;

                op++;
                op2--;
                *op = v;
                *op2 = -v;
                t++;
            }
        }
        else
        {
            while( t <= fl2 )
            {
                const double v = fvalues[ t & 1 ] *
                    pows( ( *this.*wfunc )(), Power ) / t / M_PI;

                op++;
                op2--;
                *op = v;
                *op2 = -v;
                t++;
            }
        }
    }

    template< class T >
    void generateFrac( T* op,
        CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman,
        const int opinc = 1 )
    {
        R8BASSERT( opinc != 0 );

        double f[ 2 ];
        f[ 0 ] = sin( FracDelay * M_PI );
        f[ 1 ] = -f[ 0 ];

        int t = -fl2;

        if( t + FracDelay < -Len2 )
        {
            ( *this.*wfunc )();
            *op = 0.0;
            op += opinc;
            t++;
        }

        int mt = ( FracDelay >= 1.0 - 1e-13 && FracDelay <= 1.0 + 1e-13 ?
            -1 : 0 );

        if( Power < 0.0 )
        {
            while( t < mt )
            {
                *op = f[ t & 1 ] * ( *this.*wfunc )() / ( t + FracDelay ) /
                    M_PI;

                op += opinc;
                t++;
            }

            double ut = t + FracDelay;
            *op = ( fabs( ut ) <= 1e-13 ? ( *this.*wfunc )() :
                f[ t & 1 ] * ( *this.*wfunc )() / ut / M_PI );

            mt = fl2 - 2;

            while( t < mt )
            {
                op += opinc;
                t++;
                *op = f[ t & 1 ] * ( *this.*wfunc )() / ( t + FracDelay ) /
                    M_PI;
            }

            op += opinc;
            t++;
            ut = t + FracDelay;
            *op = ( ut > Len2 ? 0.0 :
                f[ t & 1 ] * ( *this.*wfunc )() / ut / M_PI );
        }
        else
        {
            while( t < mt )
            {
                *op = f[ t & 1 ] * pows( ( *this.*wfunc )(), Power ) /
                    ( t + FracDelay ) / M_PI;

                op += opinc;
                t++;
            }

            double ut = t + FracDelay;
            *op = ( fabs( ut ) <= 1e-13 ? pows( ( *this.*wfunc )(), Power ) :
                f[ t & 1 ] * pows( ( *this.*wfunc )(), Power ) / ut / M_PI );

            mt = fl2 - 2;

            while( t < mt )
            {
                op += opinc;
                t++;
                *op = f[ t & 1 ] * pows( ( *this.*wfunc )(), Power ) /
                    ( t + FracDelay ) / M_PI;
            }

            op += opinc;
            t++;
            ut = t + FracDelay;
            *op = ( ut > Len2 ? 0.0 :
                f[ t & 1 ] * pows( ( *this.*wfunc )(), Power ) / ut / M_PI );
        }
    }

private:
    double Power; 
    CSineGen f1; 
    CSineGen f2; 
    int wn; 
    CSineGen w1; 
    CSineGen w2; 
    CSineGen w3; 

    union
    {
        struct
        {
            double KaiserBeta; 
            double KaiserDiv; 
            double KaiserLen2Frac; 
        };

        struct
        {
            double GaussianSigma; 
            double GaussianSigmaFrac; 
        };
    };

    void setWindowKaiser( const double* Params, const bool UsePower,
        const bool IsCentered )
    {
        wn = ( IsCentered ? 0 : -fl2 );

        if( Params == 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 );
        }

        KaiserDiv = besselI0( KaiserBeta );
        KaiserLen2Frac = FracDelay / Len2;
    }

    void setWindowGaussian( const double* Params, const bool UsePower,
        const bool IsCentered )
    {
        wn = ( IsCentered ? 0 : -fl2 );

        if( Params == NULL )
        {
            GaussianSigma = 1.0;
            Power = -1.0;
        }
        else
        {
            GaussianSigma = clampr( fabs( Params[ 0 ]), 1e-1, 100.0 );
            Power = ( UsePower ? fabs( Params[ 1 ]) : -1.0 );
        }

        GaussianSigma *= Len2;
        GaussianSigmaFrac = FracDelay / GaussianSigma;
    }

    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( M_PI / Len2, M_PI * 0.5 );
                w2.init( M_2PI / Len2, M_PI * 0.5 );
                w3.init( M_3PI / Len2, M_PI * 0.5 );
            }
            else
            {
                const double step1 = M_PI / Len2;
                w1.init( step1, M_PI * 0.5 - step1 * fl2 +
                    step1 * FracDelay );

                const double step2 = M_2PI / Len2;
                w2.init( step2, M_PI * 0.5 - step2 * fl2 +
                    step2 * FracDelay );

                const double step3 = M_3PI / Len2;
                w3.init( step3, M_PI * 0.5 - step3 * fl2 +
                    step3 * FracDelay );
            }

            Power = ( UsePower && Params != 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