dspop.hh

各种工程计算的库函数

/*

    DSP operations, header
    Copyright (C) 1998-2005 Jussi Laako

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

*/


#ifdef USE_INTEL_MATH
    #include <mathimf.h>
#else
    #include <math.h>
#endif
#include <float.h>

#include <Alloc.hh>

#include "dsp/dsptypes.h"
#include "dsp/DSPConfig.hh"
#include "dsp/Transform4.hh"
#include "dsp/Transform8.hh"
#include "dsp/TransformS.hh"
 
#ifdef DSP_USE_FFTW
    #include <fftw3.h>
#endif


#ifndef DSPOP_HH
    #define DSPOP_HH

    // To help compilers common subexpression elimination optimization
    #ifdef __GNUG__
        #define CONSTFUNC       __const__
    #else
        #define CONSTFUNC
    #endif
    
    // Enable/disable function inlining
    #ifndef __INTEL_COMPILER
        #define INLINE          inline
    #else
        #define INLINE
    #endif

    // function definitions for C-libraries lacking ISO C 9x standard
    #if (defined(__BCPLUSPLUS__) || defined(SOLARIS) || defined(OSX))
        #if (!defined(OSX))
            #define isgreater(x, y) (x > y)
            #define isless(x, y)    (x < y)
        #endif
        #define fmax(x, y)      ((x >= y) ? x : y)
        #define fmin(x, y)      ((x <= y) ? x : y)
        #define fmaxf(x, y)     fmax(x, y)
        #define fminf(x, y)     fmin(x, y)
        #define fma(x, y, z)    (z + x * y)
        #define fmaf(x, y, z)   fma(x, y, z)
        #define fabsf(x)        fabs(x)
        #define powf(x, y)      pow(x, y)
        #define sqrtf(x)        sqrt(x)
        #define expf(x)         exp(x)
        #define logf(x)         log(x)
        #define log10f(x)       log10(x)
        #define sinf(x)         sin(x)
        #define sinhf(x)        sinh(x)
        #define cosf(x)         cos(x)
        #define coshf(x)        cosh(x)
        #define atanf(x)        atan(x)
        #define atan2f(x, y)    atan2(x, y)
        #define acosf(x)        acos(x)
        #define acosh(x)        (log(x + sqrt(x * x - 1)))
        #define acoshf(x)       (logf(x + sqrtf(x * x - 1)))
        #define hypotf(x, y)    hypot(x, y)
    #endif

    // Microsoft VisualC++ partially conforms with ISO 9x standard
    #ifdef _MSC_VER
        #define fmax(x, y)      ((x >= y) ? x : y)
        #define fmin(x, y)      ((x <= y) ? x : y)
        #define fmaxf(x, y)     fmax(x, y)
        #define fminf(x, y)     fmin(x, y)
        #define fma(x, y, z)    (z + x * y)
        #define fmaf(x, y, z)   fma(x, y, z)
        #define acosh(x)        (log(x + sqrt(x * x - 1)))
        #define acoshf(x)       (logf(x + sqrtf(x * x - 1)))
        #define hypotf(x, y)    hypot(x, y)
        #define isgreater(x, y) (x > y)
        #define isless(x, y)    (x < y)
    #endif

    // Intel C++ partially conforms with ISO 9x standard
    #ifdef __INTEL_COMPILER
        #undef isgreater
        #undef isless
        #define isgreater(x, y) (x > y)
        #define isless(x, y)    (x < y)
    #endif

    // Maximum iterations for modified zero order Bessel function of first kind
    #define DSP_MAXBESSEL       32L

    // Filename for FFTW wisdom
    #define DSP_WISDOM_FILE     "fftw.wisdom"


    /**
        Class specialization to support automatic typecasts to cartesian/polar
        datatypes.
    */
    class clDSPAlloc : public clAlloc
    {
        public:
            clDSPAlloc () {}
            clDSPAlloc (const clDSPAlloc &CopySrc) : clAlloc(CopySrc) {}
            clDSPAlloc (long lAllocSize) : clAlloc(lAllocSize) {}
            operator stSCplx *()
                { return ((stpSCplx) vpPtr); }
            operator stDCplx *()
                { return ((stpDCplx) vpPtr); }
            operator stSPolar *()
                { return ((stpSPolar) vpPtr); }
            operator stDPolar *()
                { return ((stpDPolar) vpPtr); }
            #if defined(DSP_USE_FFTW)
            operator fftwf_complex *()
                { return ((fftwf_complex *) vpPtr); }
            operator fftw_complex *()
                { return ((fftw_complex *) vpPtr); }
            #endif
    };

    /**
        Class of DSP operations
    */
    class clDSPOp
    {
            long lPrevSrcCount;
            long lPrevDestCount;
            float fPI;
            double dPI;
            // --- Dynamically allocated arrays
            // FIR
            int iFIRDlyIdx;
            long lFIRLength;
            clDSPAlloc FIRCoeff;
            clDSPAlloc FIRBuf;
            clDSPAlloc FIRWork;
            // IIR
            float fpIIR_C[5];
            float fpIIR_X[3];
            float fpIIR_Y[2];
            double dpIIR_C[5];
            double dpIIR_X[3];
            double dpIIR_Y[2];
            // FFT (and other transforms)
            bool bFFTInitialized;
            bool bRealTransform;
            long lFFTLength;
            float fFFTScale;
            double dFFTScale;
            long *lpSBitRevWork;
            long *lpDBitRevWork;
            float *fpCosSinTable;
            double *dpCosSinTable;
            void *vpSTfrm;
            void *vpDTfrm;
            clDSPAlloc SBitRevWork;
            clDSPAlloc DBitRevWork;
            clDSPAlloc SCosSinTable;
            clDSPAlloc DCosSinTable;
            clDSPAlloc FFTBuf;
            #if defined(DSP_USE_FFTW)
            fftwf_plan fftwpSPlan;
            fftwf_plan fftwpSIPlan;
            fftw_plan fftwpDPlan;
            fftw_plan fftwpDIPlan;
            #endif
            #if defined(DSP_USE_RADIX4)
            clTransform4 Tfrm;
            #elif defined(DSP_USE_RADIX8)
            clTransform8 Tfrm;
            #else  // Split-radix
            clTransformS Tfrm;
            #endif
        protected:
            /**
                Cartesian to polar conversion.

                \f[V=\sqrt{\Re^2+\Im^2}\f]
                \f[\varphi=\arctan{\left(\frac{\Im}{\Re}\right)}\f]

                \param Magn Magnitude
                \param Phase Phase
                \param Real Real
                \param Imag Imaginary
            */
            static void Cart2Polar (float *, float *, float, float);
            /// \overload
            static void Cart2Polar (double *, double *, double, double);
            /** 
                \overload
                \param Magn Magnitude
                \param Phase Phase
                \param Cplx Cartesian
            */
            static void Cart2Polar (float *, float *, const stpSCplx);
            /// \overload
            static void Cart2Polar (double *, double *, const stpDCplx);
            /**
                \overload
                \param Polar Polar
                \param Cplx Cartesian
            */
            static void Cart2Polar (stpSPolar, const stpSCplx);
            /// \overload
            static void Cart2Polar (stpDPolar, const stpDCplx);
            /**
                \overload
                \param Coord Polar & Cartesian
            */
            static void Cart2Polar (utpSCoord);
            /// \overload
            static void Cart2Polar (utpDCoord);
            /**
                Polar to cartesian conversion.

                \f[\Re = V\cos(\varphi)\f]
                \f[\Im = V\sin(\varphi)\f]

                \note See Cart2Polar for details
            */
            static void Polar2Cart (float *, float *, float, float);
            /// \overload
            static void Polar2Cart (double *, double *, double, double);
            /// \overload
            static void Polar2Cart (stpSCplx, float, float);
            /// \overload
            static void Polar2Cart (stpDCplx, double, double);
            /// \overload
            static void Polar2Cart (stpSCplx, const stpSPolar);
            /// \overload
            static void Polar2Cart (stpDCplx, const stpDPolar);
            /// \overload
            static void Polar2Cart (utpSCoord);
            /// \overload
            static void Polar2Cart (utpDCoord);
            /**
                Complex addition

                \f[(z_{r},z_{i})=(x_{r}+y_{r},x_{i}+y_{i})\f]

                \param CplxDest Source & destination
                \param CplxSrc Source
            */
            static void CplxAdd (stpSCplx, const stpSCplx);
            /// \overload
            static void CplxAdd (stpDCplx, const stpDCplx);
            /**
                \overload
                \param CplxDest Destination
                \param CplxSrc1 Source 1
                \param CplxSrc2 Source 2
            */
            static void CplxAdd (stpSCplx, const stpSCplx, const stpSCplx);
            /// \overload
            static void CplxAdd (stpDCplx, const stpDCplx, const stpDCplx);
            /**
                Complex subtraction

                \f[(z_{r},z_{i})=(x_{r}-y_{r},x_{i}-y_{i})\f]

                \param CplxDest Source & destination
                \param CplxSrc Source
            */
            static void CplxSub (stpSCplx, const stpSCplx);
            /// \overload
            static void CplxSub (stpDCplx, const stpDCplx);
            /**
                \overload
                \param CplxDest Destination
                \param CplxSrc1 Source 1
                \param CplxSrc2 Source 2
            */
            static void CplxSub (stpSCplx, const stpSCplx, const stpSCplx);
            /// \overload
            static void CplxSub (stpDCplx, const stpDCplx, const stpDCplx);
            /**
                Complex multiply with real value

                \f[(z_r,z_i)=(x_{r}y,x_{i}y)\f]

                \param CplxDest Source & destination
                \param Src Source
            */
            static void CplxMul (stpSCplx, float);
            /// \overload
            static void CplxMul (stpDCplx, double);
            /**
                Complex multiply

                \f[\Re_{z}=\Re_{x}\Re_{y}-\Im_{x}\Im_{y}\f]
                \f[\Im_{z}=\Re_{x}\Im_{y}+\Re_{y}\Im_{x}\f]

                \param CplxDest Source & destination
                \param CplxSrc Source
            */
            static void CplxMul (stpSCplx, const stpSCplx);
            /// \overload
            static void CplxMul (stpDCplx, const stpDCplx);
            /**
                \overload
                \param CplxDest Destination
                \param CplxSrc Source 1
                \param CplxSrc Source 2
            */
            static void CplxMul (stpSCplx, const stpSCplx, const stpSCplx);
            /// \overload
            static void CplxMul (stpDCplx, const stpDCplx, const stpDCplx);
            /**
                Complex multiply with complex conjugate

                \f[\Re_{z}=\Re_{x}\Re_{y}-\Im_{x}(-\Im_{y})\f]
                \f[\Im_{z}=\Re_{x}(-\Im_{y})+\Re_{y}\Im_{x}\f]

                \param CplxDest Source & destination
                \param CplxSrc Source
            */
            static void CplxMulC (stpSCplx, const stpSCplx);
            /// \overload
            static void CplxMulC (stpDCplx, const stpDCplx);
            /**
                \overload
                \param CplxDest Destination
                \param CplxSrc1 Source 1
                \param CplxSrc2 Source 2
            */