r8butil.h

Go to the documentation of this file.

//$ nobt
//$ nocpp
 
#ifndef R8BUTIL_INCLUDED
#define R8BUTIL_INCLUDED
 
#include "r8bbase.h"
 
namespace r8b {
 
inline double convertResponseToLog( const double re, const double im )
{
    return( 4.34294481903251828 * log( re * re + im * im + 1e-100 ));
}
 
inline void updateScanStep( double& step, const double curg,
    double& prevg_log, const double prec, const double maxstep,
    const double minstep = 1e-11 )
{
    double curg_log = 4.34294481903251828 * log( curg + 1e-100 );
    curg_log += ( prevg_log - curg_log ) * 0.7;
 
    const double slope = fabs( curg_log - prevg_log );
    prevg_log = curg_log;
 
    if( slope > 0.0 )
    {
        step /= prec * slope;
        step = max( min( step, maxstep ), minstep );
    }
}
 
inline void findFIRFilterResponseMinLtoR( const double* const flt,
    const int fltlen, double& ming, double& minth, const double thend )
{
    const double maxstep = minth * 2e-3;
    double curth = minth;
    double re;
    double im;
    calcFIRFilterResponse( flt, fltlen, R8B_PI * curth, re, im );
    double prevg_log = convertResponseToLog( re, im );
    double step = 1e-11;
 
    while( true )
    {
        curth += step;
 
        if( curth > thend )
        {
            break;
        }
 
        calcFIRFilterResponse( flt, fltlen, R8B_PI * curth, re, im );
        const double curg = re * re + im * im;
 
        if( curg > ming )
        {
            ming = curg;
            minth = curth;
            break;
        }
 
        ming = curg;
        minth = curth;
 
        updateScanStep( step, curg, prevg_log, 0.31, maxstep );
    }
}
 
inline void findFIRFilterResponseMaxLtoR( const double* const flt,
    const int fltlen, double& maxg, double& maxth, const double thend )
{
    const double maxstep = maxth * 1e-4;
    double premaxth = maxth;
    double premaxg = maxg;
    double postmaxth = maxth;
    double postmaxg = maxg;
 
    double prevth = maxth;
    double prevg = maxg;
    double curth = maxth;
    double re;
    double im;
    calcFIRFilterResponse( flt, fltlen, R8B_PI * curth, re, im );
    double prevg_log = convertResponseToLog( re, im );
    double step = 1e-11;
 
    bool WasPeak = false;
    int AfterPeakCount = 0;
 
    while( true )
    {
        curth += step;
 
        if( curth > thend )
        {
            break;
        }
 
        calcFIRFilterResponse( flt, fltlen, R8B_PI * curth, re, im );
        const double curg = re * re + im * im;
 
        if( curg > maxg )
        {
            premaxth = prevth;
            premaxg = prevg;
            maxg = curg;
            maxth = curth;
            WasPeak = true;
            AfterPeakCount = 0;
        }
        else
        if( WasPeak )
        {
            if( AfterPeakCount == 0 )
            {
                postmaxth = curth;
                postmaxg = curg;
            }
 
            if( AfterPeakCount == 5 )
            {
                // Perform 2 approximate binary searches.
 
                int k;
 
                for( k = 0; k < 2; k++ )
                {
                    double l = ( k == 0 ? premaxth : maxth );
                    double curgl = ( k == 0 ? premaxg : maxg );
                    double r = ( k == 0 ? maxth : postmaxth );
                    double curgr = ( k == 0 ? maxg : postmaxg );
 
                    while( true )
                    {
                        const double c = ( l + r ) * 0.5;
                        calcFIRFilterResponse( flt, fltlen, R8B_PI * c,
                            re, im );
 
                        const double curg = re * re + im * im;
 
                        if( curgl > curgr )
                        {
                            r = c;
                            curgr = curg;
                        }
                        else
                        {
                            l = c;
                            curgl = curg;
                        }
 
                        if( r - l < 1e-11 )
                        {
                            if( curgl > curgr )
                            {
                                maxth = l;
                                maxg = curgl;
                            }
                            else
                            {
                                maxth = r;
                                maxg = curgr;
                            }
 
                            break;
                        }
                    }
                }
 
                break;
            }
 
            AfterPeakCount++;
        }
 
        prevth = curth;
        prevg = curg;
 
        updateScanStep( step, curg, prevg_log, 1.0, maxstep );
    }
}
 
inline void findFIRFilterResponseLevelRtoL( const double* const flt,
    const int fltlen, const double maxg, double& th, const double thend )
{
    // Perform exact binary search.
 
    double l = thend;
    double r = th;
 
    while( true )
    {
        const double c = ( l + r ) * 0.5;
 
        if( r - l < 1e-14 )
        {
            th = c;
            break;
        }
 
        double re;
        double im;
        calcFIRFilterResponse( flt, fltlen, R8B_PI * c, re, im );
        const double curg = re * re + im * im;
 
        if( curg > maxg )
        {
            l = c;
        }
        else
        {
            r = c;
        }
    }
}
 
} // namespace r8b
 
#endif // R8BUTIL_INCLUDED