CDSPFIRFilter.h

Go to the documentation of this file.

//$ nobt
//$ nocpp

 
#ifndef R8B_CDSPFIRFILTER_INCLUDED
#define R8B_CDSPFIRFILTER_INCLUDED
 
#include "CDSPSincFilterGen.h"
#include "CDSPRealFFT.h"
 
namespace r8b {

 
enum EDSPFilterPhaseResponse
{
    fprLinearPhase = 0, 
    fprMinPhase 
};

 
class CDSPFIRFilter : public R8B_BASECLASS,
    private CSinglyLinkedListItem< CDSPFIRFilter >
{
    R8BNOCTOR( CDSPFIRFilter )
 
    friend class CSinglyLinkedListItem< CDSPFIRFilter >;
    friend class CDSPFIRFilterCache;
 
public:
 
    static double getLPMinTransBand()
    {
        return( 0.5 );
    }

 
    static double getLPMaxTransBand()
    {
        return( 45.0 );
    }

 
    static double getLPMinAtten()
    {
        return( 49.0 );
    }

 
    static double getLPMaxAtten()
    {
        return( 218.0 );
    }

 
    bool isZeroPhase() const
    {
        return( IsZeroPhase );
    }

 
    int getLatency() const
    {
        return( Latency );
    }

 
    double getLatencyFrac() const
    {
        return( LatencyFrac );
    }

 
    int getKernelLen() const
    {
        return( KernelLen );
    }

 
    int getBlockLenBits() const
    {
        return( BlockLenBits );
    }

 
    const realfft_t* getKernelBlock() const
    {
        return( KernelBlockFFT );
    }

 
    realfft_t calcDownShiftNyquist( const realfft_t re,
        const realfft_t im ) const
    {
        return( dsnRe * re - dsnIm * im );
    }

 
    void unref();
 
private:
    double ReqNormFreq; 
    double ReqTransBand; 
    double ReqAtten; 
    EDSPFilterPhaseResponse ReqPhase; 
    double ReqGain; 
    int ReqDownShift; 
    int RefCount; 
    bool IsZeroPhase; 
    int Latency; 
    double LatencyFrac; 
    int KernelLen; 
    int BlockLenBits; 
    CFixedBuffer< double > KernelBlock; 
    realfft_t* KernelBlockFFT; 
    realfft_t dsnRe; 
    realfft_t dsnIm; 
 
    CDSPFIRFilter()
        : RefCount( 1 )
    {
    }

 
    void buildLPFilter( const double* const ExtAttenCorrs )
    {
        const double tb = ReqTransBand * 0.01;
        double pwr;
        double fo1;
        double hl;
        double atten = -ReqAtten;
 
        if( tb >= 0.25 )
        {
            if( ReqAtten >= 117.0 )
            {
                atten -= 1.60;
            }
            else
            if( ReqAtten >= 60.0 )
            {
                atten -= 1.91;
            }
            else
            {
                atten -= 2.25;
            }
        }
        else
        if( tb >= 0.10 )
        {
            if( ReqAtten >= 117.0 )
            {
                atten -= 0.69;
            }
            else
            if( ReqAtten >= 60.0 )
            {
                atten -= 0.73;
            }
            else
            {
                atten -= 1.13;
            }
        }
        else
        {
            if( ReqAtten >= 117.0 )
            {
                atten -= 0.21;
            }
            else
            if( ReqAtten >= 60.0 )
            {
                atten -= 0.25;
            }
            else
            {
                atten -= 0.36;
            }
        }
 
        static const int AttenCorrCount = 264;
        static const double AttenCorrMin = 49.0;
        static const double AttenCorrDiff = 176.25;
        int AttenCorr = (int) floor(( -atten - AttenCorrMin ) *
            AttenCorrCount / AttenCorrDiff + 0.5 );
 
        AttenCorr = min( AttenCorrCount, max( 0, AttenCorr ));
 
        if( ExtAttenCorrs != R8B_NULL )
        {
            atten -= ExtAttenCorrs[ AttenCorr ];
        }
        else
        if( tb >= 0.25 )
        {
            static const double AttenCorrScale = 101.0;
            static const signed char AttenCorrs[] = {
                -127, -127, -125, -125, -122, -119, -115, -110, -104, -97,
                -91, -82, -75, -24, -16, -6, 4, 14, 24, 29, 30, 32, 37, 44,
                51, 57, 63, 67, 65, 50, 53, 56, 58, 60, 63, 64, 66, 68, 74,
                77, 78, 78, 78, 79, 79, 60, 60, 60, 61, 59, 52, 47, 41, 36,
                30, 24, 17, 9, 0, -8, -10, -11, -14, -13, -18, -25, -31, -38,
                -44, -50, -57, -63, -68, -74, -81, -89, -96, -101, -104, -107,
                -109, -110, -86, -84, -85, -82, -80, -77, -73, -67, -62, -55,
                -48, -42, -35, -30, -20, -11, -2, 5, 6, 6, 7, 11, 16, 21, 26,
                34, 41, 46, 49, 52, 55, 56, 48, 49, 51, 51, 52, 52, 52, 52,
                52, 51, 51, 50, 47, 47, 50, 48, 46, 42, 38, 35, 31, 27, 24,
                20, 16, 12, 11, 12, 10, 8, 4, -1, -6, -11, -16, -19, -17, -21,
                -24, -27, -32, -34, -37, -38, -40, -41, -40, -40, -42, -41,
                -44, -45, -43, -41, -34, -31, -28, -24, -21, -18, -14, -10,
                -5, -1, 2, 5, 8, 7, 4, 3, 2, 2, 4, 6, 8, 9, 9, 10, 10, 10, 10,
                9, 8, 9, 11, 14, 13, 12, 11, 10, 8, 7, 6, 5, 3, 2, 2, -1, -1,
                -3, -3, -4, -4, -5, -4, -6, -7, -9, -5, -1, -1, 0, 1, 0, -2,
                -3, -4, -5, -5, -8, -13, -13, -13, -12, -13, -12, -11, -11,
                -9, -8, -7, -5, -3, -1, 2, 4, 6, 9, 10, 11, 14, 18, 21, 24,
                27, 30, 34, 37, 37, 39, 40 };
 
            atten -= AttenCorrs[ AttenCorr ] / AttenCorrScale;
        }
        else
        if( tb >= 0.10 )
        {
            static const double AttenCorrScale = 210.0;
            static const signed char AttenCorrs[] = {
                -113, -118, -122, -125, -126, -97, -95, -92, -92, -89, -82,
                -75, -69, -48, -42, -36, -30, -22, -14, -5, -2, 1, 6, 13, 22,
                28, 35, 41, 48, 55, 56, 56, 61, 65, 71, 77, 81, 83, 85, 85,
                74, 74, 73, 72, 71, 70, 68, 64, 59, 56, 49, 52, 46, 42, 36,
                32, 26, 20, 13, 7, -2, -6, -10, -15, -20, -27, -33, -38, -44,
                -43, -48, -53, -57, -63, -69, -73, -75, -79, -81, -74, -76,
                -77, -77, -78, -81, -80, -80, -78, -76, -65, -62, -59, -56,
                -51, -48, -44, -38, -33, -25, -19, -13, -5, -1, 2, 7, 13, 17,
                21, 25, 30, 35, 40, 45, 50, 53, 56, 57, 55, 58, 59, 62, 64,
                67, 67, 68, 68, 62, 61, 61, 59, 59, 57, 57, 55, 52, 48, 42,
                38, 35, 31, 26, 20, 15, 13, 10, 7, 3, -2, -8, -13, -17, -23,
                -28, -34, -37, -40, -41, -45, -48, -50, -53, -57, -59, -62,
                -63, -63, -57, -57, -56, -56, -54, -54, -53, -49, -48, -41,
                -38, -33, -31, -26, -23, -18, -12, -9, -7, -7, -3, 0, 5, 9,
                14, 16, 20, 22, 21, 23, 25, 27, 28, 29, 34, 33, 35, 33, 31,
                30, 29, 29, 26, 26, 25, 24, 20, 19, 15, 10, 8, 4, 1, -2, -6,
                -10, -16, -19, -23, -26, -27, -30, -34, -39, -43, -47, -51,
                -52, -54, -56, -58, -59, -62, -63, -66, -65, -65, -64, -59,
                -57, -54, -52, -48, -44, -42, -37, -32, -22, -17, -10, -3, 5,
                13, 22, 30, 40, 50, 60, 72 };
 
            atten -= AttenCorrs[ AttenCorr ] / AttenCorrScale;
        }
        else
        {
            static const double AttenCorrScale = 196.0;
            static const signed char AttenCorrs[] = {
                -15, -17, -20, -20, -20, -21, -20, -16, -17, -18, -17, -13,
                -12, -11, -9, -7, -5, -4, -1, 1, 3, 4, 5, 6, 7, 9, 9, 10, 10,
                10, 11, 11, 11, 12, 12, 12, 10, 11, 10, 10, 8, 10, 11, 10, 11,
                11, 13, 14, 15, 19, 27, 26, 23, 18, 14, 8, 4, -2, -6, -12,
                -17, -23, -28, -33, -37, -42, -46, -49, -53, -57, -60, -61,
                -64, -65, -67, -66, -66, -66, -65, -64, -61, -59, -56, -52,
                -48, -42, -38, -31, -27, -19, -13, -7, -1, 8, 14, 22, 29, 37,
                45, 52, 59, 66, 73, 80, 86, 91, 96, 100, 104, 108, 111, 114,
                115, 117, 118, 120, 120, 118, 117, 114, 113, 111, 107, 103,
                99, 95, 89, 84, 78, 72, 66, 60, 52, 44, 37, 30, 21, 14, 6, -3,
                -11, -18, -26, -34, -43, -51, -58, -65, -73, -78, -85, -90,
                -97, -102, -107, -113, -115, -118, -121, -125, -125, -126,
                -126, -126, -125, -124, -121, -119, -115, -111, -109, -101,
                -102, -95, -88, -81, -73, -67, -63, -54, -47, -40, -33, -26,
                -18, -11, -5, 2, 8, 14, 19, 25, 31, 36, 37, 43, 47, 49, 51,
                52, 57, 57, 56, 57, 58, 58, 58, 57, 56, 52, 52, 50, 48, 44,
                41, 39, 37, 33, 31, 26, 24, 21, 18, 14, 11, 8, 4, 2, -2, -5,
                -7, -9, -11, -13, -15, -16, -18, -19, -20, -23, -24, -24, -25,
                -27, -26, -27, -29, -30, -31, -32, -35, -36, -39, -40, -44,
                -46, -51, -54, -59, -63, -69, -76, -83, -91, -98 };
 
            atten -= AttenCorrs[ AttenCorr ] / AttenCorrScale;
        }
 
        pwr = 7.43932822146293e-8 * sqr( atten ) + 0.000102747434588003 *
            cos( 0.00785021930010397 * atten ) * cos( 0.633854318781239 +
            0.103208573657699 * atten ) - 0.00798132247867036 -
            0.000903555213543865 * atten - 0.0969365532127236 * exp(
            0.0779275237937911 * atten ) - 1.37304948662012e-5 * atten * cos(
            0.00785021930010397 * atten );
 
        if( pwr <= 0.067665322581 )
        {
            if( tb >= 0.25 )
            {
                hl = 2.6778150875894 / tb + 300.547590563091 * atan( atan(
                    2.68959772209918 * pwr )) / ( 5.5099277187035 * tb - tb *
                    tanh( cos( asinh( atten ))));
 
                fo1 = 0.987205355829873 * tb + 1.00011788929851 * atan2(
                    -0.321432067051302 - 6.19131357321578 * sqrt( pwr ),
                    hl + -1.14861472207245 / ( hl - 14.1821147585957 ) + pow(
                    0.9521145021664, pow( atan2( 1.12018764830637, tb ),
                    2.10988901686912 * hl - 20.9691278378345 )));
            }
            else
            if( tb >= 0.10 )
            {
                hl = ( 1.56688617018066 + 142.064321294568 * pwr +
                    0.00419441117131136 * cos( 243.633511747297 * pwr ) -
                    0.022953443903576 * atten - 0.026629568860284 * cos(
                    127.715550622571 * pwr )) / tb;
 
                fo1 = 0.982299356642411 * tb + 0.999441744774215 * asinh((
                    -0.361783054039583 - 5.80540593623676 * sqrt( pwr )) /
                    hl );
            }
            else
            {
                hl = ( 2.45739657014937 + 269.183679500541 * pwr * cos(
                    5.73225668178813 + atan2( cosh( 0.988861169868941 -
                    17.2201556280744 * pwr ), 1.08340138240431 * pwr ))) / tb;
 
                fo1 = 2.291956939 * tb + 0.01942450693 * sqr( tb ) * hl -
                    4.67538973161837 * pwr * tb - 1.668433124 * tb *
                    pow( pwr, pwr );
            }
        }
        else
        {
            if( tb >= 0.25 )
            {
                hl = ( 1.50258368698213 + 158.556968859477 * asinh( pwr ) *
                    tanh( 57.9466246871383 * tanh( pwr )) -
                    0.0105440479814834 * atten ) / tb;
 
                fo1 = 0.994024401639321 * tb + ( -0.236282717577215 -
                    6.8724924545387 * sqrt( sin( pwr ))) / hl;
            }
            else
            if( tb >= 0.10 )
            {
                hl = ( 1.50277377248945 + 158.222625721046 * asinh( pwr ) *
                    tanh( 1.02875299001715 + 42.072277322604 * pwr ) -
                    0.0108380943845632 * atten ) / tb;
 
                fo1 = 0.992539376734551 * tb + ( -0.251747813037178 -
                    6.74159892452584 * sqrt( tanh( tanh( tan( pwr ))))) / hl;
            }
            else
            {
                hl = ( 1.15990238966306 * pwr - 5.02124037125213 * sqr(
                    pwr ) - 0.158676856669827 * atten * cos( 1.1609073390614 *
                    pwr - 6.33932586197475 * pwr * sqr( pwr ))) / tb;
 
                fo1 = 0.867344453126885 * tb + 0.052693817907757 * tb * log(
                    pwr ) + 0.0895511178735932 * tb * atan( 59.7538527741309 *
                    pwr ) - 0.0745653568081453 * pwr * tb;
            }
        }
 
        double WinParams[ 2 ];
        WinParams[ 0 ] = 125.0;
        WinParams[ 1 ] = pwr;
 
        CDSPSincFilterGen sinc;
        sinc.Len2 = 0.25 * hl / ReqNormFreq;
        sinc.Freq1 = 0.0;
        sinc.Freq2 = R8B_PI * ( 1.0 - fo1 ) * ReqNormFreq;
        sinc.initBand( CDSPSincFilterGen :: wftKaiser, WinParams, true );
 
        KernelLen = sinc.KernelLen;
        BlockLenBits = getBitOccupancy( KernelLen - 1 ) + R8B_EXTFFT;
        const int BlockLen = 1 << BlockLenBits;
 
        KernelBlock.alloc( BlockLen * 2 );
        sinc.generateBand( &KernelBlock[ 0 ],
            &CDSPSincFilterGen :: calcWindowKaiser );
 
        if( ReqPhase == fprLinearPhase )
        {
            IsZeroPhase = true;
            Latency = sinc.fl2;
            LatencyFrac = 0.0;
        }
        else
        {
            IsZeroPhase = false;
            double DCGroupDelay;
 
            calcMinPhaseTransform( &KernelBlock[ 0 ], KernelLen, 16, false,
                &DCGroupDelay );
 
            Latency = (int) DCGroupDelay;
            LatencyFrac = DCGroupDelay - Latency;
        }
 
        CDSPRealFFTKeeper ffto( BlockLenBits + 1 );
 
        dsnRe = 0;
        dsnIm = 0;
 
        if( IsZeroPhase )
        {
            // Calculate DC gain.
 
            double s = 0.0;
            int i;
 
            for( i = 0; i < KernelLen; i++ )
            {
                s += KernelBlock[ i ];
            }
 
            s = ffto -> getInvMulConst() * ReqGain / s;
 
            // Time-shift the filter so that zero-phase response is produced.
            // Simultaneously multiply by "s".
 
            for( i = 0; i <= sinc.fl2; i++ )
            {
                KernelBlock[ i ] = KernelBlock[ sinc.fl2 + i ] * s;
            }
 
            for( i = 1; i <= sinc.fl2; i++ )
            {
                KernelBlock[ BlockLen * 2 - i ] = KernelBlock[ i ];
            }
 
            memset( &KernelBlock[ sinc.fl2 + 1 ], 0,
                (size_t) ( BlockLen * 2 - KernelLen ) *
                sizeof( KernelBlock[ 0 ]));
 
            KernelBlockFFT = ffto -> forward( KernelBlock );
 
            if( ReqDownShift > 0 )
            {
                KernelBlockFFT = ffto -> reorderForward( KernelBlockFFT,
                    ffto -> getWorkBuf() );
 
                const int z = ( 2 << BlockLenBits ) >> ReqDownShift;
                dsnRe = KernelBlockFFT[ z ];
                dsnIm = KernelBlockFFT[ z + 1 ];
 
                CDSPRealFFTKeeper ffto2( BlockLenBits + 1 - ReqDownShift );
                CFixedBuffer< double > NewBlock( ffto2 -> getLen() *
                    (int) sizeof( realfft_t ) / (int) sizeof( double ));
 
                KernelBlockFFT = ffto2 -> reorderInverse( KernelBlockFFT,
                    NewBlock );
 
                ffto2 -> convertToZP( KernelBlockFFT );
 
                KernelBlock.moveFrom( NewBlock );
            }
            else
            {
                ffto -> convertToZP( KernelBlockFFT );
            }
        }
        else
        {
            normalizeFIRFilter( &KernelBlock[ 0 ], KernelLen,
                ffto -> getInvMulConst() * ReqGain );
 
            memset( &KernelBlock[ KernelLen ], 0,
                (size_t) ( BlockLen * 2 - KernelLen ) *
                sizeof( KernelBlock[ 0 ]));
 
            KernelBlockFFT = ffto -> forward( KernelBlock );
 
            if( ReqDownShift > 0 )
            {
                KernelBlockFFT = ffto -> reorderForward( KernelBlockFFT,
                    ffto -> getWorkBuf() );
 
                const int z = ( 2 << BlockLenBits ) >> ReqDownShift;
                dsnRe = KernelBlockFFT[ z ];
                dsnIm = KernelBlockFFT[ z + 1 ];
 
                CDSPRealFFTKeeper ffto2( BlockLenBits + 1 - ReqDownShift );
                CFixedBuffer< double > NewBlock( ffto2 -> getLen() *
                    (int) sizeof( realfft_t ) / (int) sizeof( double ));
 
                KernelBlockFFT = ffto2 -> reorderInverse( KernelBlockFFT,
                    NewBlock );
 
                KernelBlock.moveFrom( NewBlock );
            }
        }
 
        R8BCONSOLE( "CDSPFIRFilter: flt_len=%i latency=%i nfreq=%.4f "
            "tb=%.1f att=%.1f gain=%.3f dsh=%i\n", KernelLen, Latency,
            ReqNormFreq, ReqTransBand, ReqAtten, ReqGain, ReqDownShift );
    }
};

 
class CDSPFIRFilterCache : public R8B_BASECLASS
{
    R8BNOCTOR( CDSPFIRFilterCache )
 
    friend class CDSPFIRFilter;
 
public:
 
    static int getObjCount()
    {
        R8BSYNC( getStateSync() );
 
        return( getObjCountStatic() );
    }

 
    static CDSPFIRFilter& getLPFilter( const double ReqNormFreq,
        const double ReqTransBand, const double ReqAtten,
        const EDSPFilterPhaseResponse ReqPhase, const double ReqGain,
        const int ReqDownFactor, const double* const AttenCorrs = R8B_NULL )
    {
        R8BASSERT( ReqNormFreq > 0.0 && ReqNormFreq <= 1.0 );
        R8BASSERT( ReqTransBand >= CDSPFIRFilter :: getLPMinTransBand() );
        R8BASSERT( ReqTransBand <= CDSPFIRFilter :: getLPMaxTransBand() );
        R8BASSERT( ReqAtten >= CDSPFIRFilter :: getLPMinAtten() );
        R8BASSERT( ReqAtten <= CDSPFIRFilter :: getLPMaxAtten() );
        R8BASSERT( ReqGain > 0.0 );
        R8BASSERT( ReqDownFactor > 0 );
 
        R8B_EXITDTOR static CPtrKeeper< CDSPFIRFilter > Objects; // The chain
            // of cached objects.
 
        int ReqDownShift = getBitOccupancy( ReqDownFactor ) - 1;
 
        if(( 1 << ReqDownShift ) != ReqDownFactor )
        {
            ReqDownShift = 0;
        }
 
        R8BSYNC( getStateSync() );
 
        int& ObjCount = getObjCountStatic();
        CDSPFIRFilter* PrevObj = R8B_NULL;
        CDSPFIRFilter* CurObj = Objects;
 
        while( CurObj != R8B_NULL )
        {
            if( CurObj -> ReqNormFreq == ReqNormFreq &&
                CurObj -> ReqTransBand == ReqTransBand &&
                CurObj -> ReqGain == ReqGain &&
                CurObj -> ReqDownShift == ReqDownShift &&
                CurObj -> ReqAtten == ReqAtten &&
                CurObj -> ReqPhase == ReqPhase )
            {
                break;
            }
 
            if( CurObj -> Next == R8B_NULL &&
                ObjCount >= R8B_FILTER_CACHE_MAX )
            {
                if( CurObj -> RefCount == 0 )
                {
                    // Delete the last filter which is not used.
 
                    PrevObj -> Next = R8B_NULL;
                    delete CurObj;
                    ObjCount--;
                }
                else
                {
                    // Move the last filter to the top of the list since it
                    // seems to be in use for a long time.
 
                    PrevObj -> Next = R8B_NULL;
                    CurObj -> Next = Objects.unkeep();
                    Objects = CurObj;
                }
 
                CurObj = R8B_NULL;
                break;
            }
 
            PrevObj = CurObj;
            CurObj = CurObj -> Next;
        }
 
        if( CurObj != R8B_NULL )
        {
            CurObj -> RefCount++;
 
            if( PrevObj == R8B_NULL )
            {
                return( *CurObj );
            }
 
            // Remove the filter from the list temporarily.
 
            PrevObj -> Next = CurObj -> Next;
        }
        else
        {
            // Create a new filter object (with RefCount == 1) and build the
            // filter kernel.
 
            CPtrKeeper< CDSPFIRFilter > f( new CDSPFIRFilter() );
 
            f -> ReqNormFreq = ReqNormFreq;
            f -> ReqTransBand = ReqTransBand;
            f -> ReqAtten = ReqAtten;
            f -> ReqPhase = ReqPhase;
            f -> ReqGain = ReqGain;
            f -> ReqDownShift = ReqDownShift;
 
            f -> buildLPFilter( AttenCorrs );
 
            CurObj = f.unkeep();
            ObjCount++;
        }
 
        // Insert the filter at the start of the list.
 
        CurObj -> Next = Objects.unkeep();
        Objects = CurObj;
 
        return( *CurObj );
    }
 
protected:
 
    static CSyncObject& getStateSync()
    {
        R8B_EXITDTOR static CSyncObject StateSync;
 
        return( StateSync );
    }

 
    static int& getObjCountStatic()
    {
        R8B_EXITDTOR static int ObjCount = 0; // The number of objects
            // currently preset in the cache.
 
        return( ObjCount );
    }
};
 
// ---------------------------------------------------------------------------
// CDSPFIRFilter PUBLIC
// ---------------------------------------------------------------------------
 
inline void CDSPFIRFilter :: unref()
{
    R8BSYNC( CDSPFIRFilterCache :: getStateSync() );
 
    RefCount--;
}
 
// ---------------------------------------------------------------------------
 
} // namespace r8b
 
#endif // R8B_CDSPFIRFILTER_INCLUDED