dspop.cc

各种工程计算的库函数

    dRealExp = exp(spCplxSrc->R);
    #ifndef _GNU_SOURCE
    spCplxDest->R = dRealExp * cos(spCplxSrc->I);
    spCplxDest->I = dRealExp * sin(spCplxSrc->I);
    #else
    sincos(spCplxSrc->I, &spCplxDest->I, &spCplxDest->R);
    spCplxDest->R *= dRealExp;
    spCplxDest->I *= dRealExp;
    #endif
}


INLINE void clDSPOp::CplxLog (stpSCplx spCplxDest, const stpSCplx spCplxSrc)
{
    Cart2Polar(&spCplxDest->R, &spCplxDest->I, spCplxSrc->R, spCplxSrc->I);
    spCplxDest->R = logf(spCplxDest->R);
}


INLINE void clDSPOp::CplxLog (stpDCplx spCplxDest, const stpDCplx spCplxSrc)
{
    Cart2Polar(&spCplxDest->R, &spCplxDest->I, spCplxSrc->R, spCplxSrc->I);
    spCplxDest->R = log(spCplxDest->R);
}


INLINE void clDSPOp::CplxLog10 (stpSCplx spCplxDest, const stpSCplx spCplxSrc)
{
    Cart2Polar(&spCplxDest->R, &spCplxDest->I, spCplxSrc->R, spCplxSrc->I);
    spCplxDest->R = log10f(spCplxDest->R);
}


INLINE void clDSPOp::CplxLog10 (stpDCplx spCplxDest, const stpDCplx spCplxSrc)
{
    Cart2Polar(&spCplxDest->R, &spCplxDest->I, spCplxSrc->R, spCplxSrc->I);
    spCplxDest->R = log10(spCplxDest->R);
}


INLINE void clDSPOp::CplxPow (stpSCplx spCplxDest, const stpSCplx spCplxSrc,
    const stpSCplx spCplxExp)
{
    stSCplx sCplxTemp;

    CplxLog(&sCplxTemp, spCplxSrc);
    CplxMul(&sCplxTemp, spCplxExp);
    CplxExp(spCplxDest, &sCplxTemp);
}


INLINE void clDSPOp::CplxPow (stpDCplx spCplxDest, const stpDCplx spCplxSrc,
    const stpDCplx spCplxExp)
{
    stDCplx sCplxTemp;

    CplxLog(&sCplxTemp, spCplxSrc);
    CplxMul(&sCplxTemp, spCplxExp);
    CplxExp(spCplxDest, &sCplxTemp);
}


INLINE void clDSPOp::CplxRoot (stpSCplx spCplxDest, const stpSCplx spCplxSrc,
    const stpSCplx spCplxRoot)
{
    stSCplx sCplxExp;

    sCplxExp.R = 1.0F;
    sCplxExp.I = 0.0F;
    CplxDiv(&sCplxExp, spCplxRoot);
    CplxPow(spCplxDest, spCplxSrc, &sCplxExp);
}


INLINE void clDSPOp::CplxRoot (stpDCplx spCplxDest, const stpDCplx spCplxSrc,
    const stpDCplx spCplxRoot)
{
    stDCplx sCplxExp;

    sCplxExp.R = 1.0;
    sCplxExp.I = 0.0;
    CplxDiv(&sCplxExp, spCplxRoot);
    CplxPow(spCplxDest, spCplxSrc, &sCplxExp);
}


/*INLINE void clDSPOp::CplxConj (stpSCplx spCplx)
{
    spCplx->I = -(spCplx->I);
}


INLINE void clDSPOp::CplxConj (stpDCplx spCplx)
{
    spCplx->I = -(spCplx->I);
}*/


INLINE void clDSPOp::CplxConj (stpSCplx spCplxDest, const stpSCplx spCplxSrc)
{
    spCplxDest->R = spCplxSrc->R;
    spCplxDest->I = -(spCplxSrc->I);
}


INLINE void clDSPOp::CplxConj (stpDCplx spCplxDest, const stpDCplx spCplxSrc)
{
    spCplxDest->R = spCplxSrc->R;
    spCplxDest->I = -(spCplxSrc->I);
}


INLINE double clDSPOp::Multiple (long lValue)
{
    long lLoopCntr;
    double dMult = 1.0;

    for (lLoopCntr = 1L; lLoopCntr <= lValue; lLoopCntr++)
    {
        dMult *= (double) lLoopCntr;
    }
    return dMult;
}


INLINE float clDSPOp::ModZeroBessel (float fValue)
{
    long lLoopCntr;
    float fMZBessel = 0.0F;
    float fHalfValue;

    fHalfValue = fValue / 2.0F;
    for (lLoopCntr = 0L; lLoopCntr <= DSP_MAXBESSEL; lLoopCntr++)
    {
        fMZBessel += (float) pow(
            pow(fHalfValue, lLoopCntr) / Multiple(lLoopCntr), 2.0);
    }
    return fMZBessel;
}


INLINE double clDSPOp::ModZeroBessel (double dValue)
{
    long lLoopCntr;
    double dMZBessel = 0.0;
    double dHalfValue;

    dHalfValue = dValue / 2.0;
    for (lLoopCntr = 0L; lLoopCntr <= DSP_MAXBESSEL; lLoopCntr++)
    {
        dMZBessel += pow(
            pow(dHalfValue, lLoopCntr) / Multiple(lLoopCntr), 2.0);
    }
    return dMZBessel;
}


float clDSPOp::ChebyshevPolynom (float fOrder, float fValue)
{
    if (fabsf(fValue) <= 1.0F)
    {
        return cosf(fOrder * acosf(fValue));
    }
    else
    {
        return coshf(fOrder * acoshf(fValue));
    }
}


double clDSPOp::ChebyshevPolynom (double dOrder, double dValue)
{
    if (fabs(dValue) <= 1.0)
    {
        return cos(dOrder * acos(dValue));
    }
    else
    {
        return cosh(dOrder * acosh(dValue));
    }
}


clDSPOp::clDSPOp()
{
    #ifdef DSP_X86
    bHave3DNow = (dsp_x86_have_e3dnow()) ? true : false;
    bHaveSSE = (dsp_x86_have_sse2()) ? true : false;
    #ifdef __INTEL_COMPILER
    if (bHaveSSE)
    {
        _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
        _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
    }
    #endif
    #endif
    lPrevSrcCount = 0;
    lPrevDestCount = 0;
    // This is to get maximum precision for PI
    #ifdef _ISOC9X_SOURCE
    fPI = (float) acosl(-1.0);
    dPI = (double) acosl(-1.0);
    #else
    fPI = (float) acos(-1.0);
    dPI = acos(-1.0);
    #endif
    bFFTInitialized = false;
    lpSBitRevWork = NULL;
    lpDBitRevWork = NULL;
    fpCosSinTable = NULL;
    dpCosSinTable = NULL;
}


clDSPOp::~clDSPOp()
{
    if (bFFTInitialized)
    {
        FFTUninitialize();
    }
}


void clDSPOp::Add (float *fpDest, float fSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsAddC_32f_I(fSrc, fpDest, lCount);
    #else
    long lLoopCntr;

    #ifdef DSP_X86
    if (bHave3DNow)
    {
        dsp_x86_3dnow_addf(fpDest, fSrc, lCount);
    }
    else if (bHaveSSE)
    {
        dsp_x86_sse_addf(fpDest, fSrc, lCount);
    }
    else
    #endif
    {
        for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
        {
            fpDest[lLoopCntr] += fSrc;
        }
    }
    #endif
}


void clDSPOp::Add (double *dpDest, double dSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsAddC_64f_I(dSrc, dpDest, lCount);
    #else
    long lLoopCntr;

    #ifdef DSP_X86
    if (bHaveSSE)
    {
        dsp_x86_sse_add(dpDest, dSrc, lCount);
    }
    else
    #endif
    {
        for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
        {
            dpDest[lLoopCntr] += dSrc;
        }
    }
    #endif
}


void clDSPOp::Add (stpSCplx spCplxDest, stSCplx sCplxSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsAddC_32fc_I(*((Ipp32fc *) &sCplxSrc), (Ipp32fc *) spCplxDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        CplxAdd(&spCplxDest[lLoopCntr], &sCplxSrc);
    }
    #endif
}


void clDSPOp::Add (stpDCplx spCplxDest, stDCplx sCplxSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsAddC_64fc_I(*((Ipp64fc *) &sCplxSrc), (Ipp64fc *) spCplxDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        CplxAdd(&spCplxDest[lLoopCntr], &sCplxSrc);
    }
    #endif
}


void clDSPOp::Add (float *fpDest, const float *fpSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsAdd_32f_I(fpSrc, fpDest, lCount);
    #else
    long lLoopCntr;

    #ifdef DSP_X86
    if (bHave3DNow)
    {
        dsp_x86_3dnow_add2f(fpDest, fpSrc, lCount);
    }
    else if (bHaveSSE)
    {
        dsp_x86_sse_add2f(fpDest, fpSrc, lCount);
    }
    else
    #endif
    {
        for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
        {
            fpDest[lLoopCntr] += fpSrc[lLoopCntr];
        }
    }
    #endif
}


void clDSPOp::Add (double *dpDest, const double *dpSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsAdd_64f_I(dpSrc, dpDest, lCount);
    #else
    long lLoopCntr;

    #ifdef DSP_X86
    if (bHaveSSE)
    {
        dsp_x86_sse_add2(dpDest, dpSrc, lCount);
    }
    else
    #endif
    {
        for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
        {
            dpDest[lLoopCntr] += dpSrc[lLoopCntr];
        }
    }
    #endif
}


void clDSPOp::Add (stpSCplx spCplxDest, const stpSCplx spCplxSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsAdd_32fc_I((Ipp32fc *) spCplxSrc, (Ipp32fc *) spCplxDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        CplxAdd(&spCplxDest[lLoopCntr], &spCplxSrc[lLoopCntr]);
    }
    #endif
}


void clDSPOp::Add (stpDCplx spCplxDest, const stpDCplx spCplxSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsAdd_64fc_I((Ipp64fc *) spCplxSrc, (Ipp64fc *) spCplxDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        CplxAdd(&spCplxDest[lLoopCntr], &spCplxSrc[lLoopCntr]);
    }
    #endif
}


void clDSPOp::Add (float *fpDest, const float *fpSrc1, const float *fpSrc2, 
    long lCount)
{
    #ifdef DSP_IPP
    ippsAdd_32f(fpSrc1, fpSrc2, fpDest, lCount);
    #else
    long lLoopCntr;

    #ifdef DSP_X86
    if (bHave3DNow)
    {
        dsp_x86_3dnow_add3f(fpDest, fpSrc1, fpSrc2, lCount);
    }
    else if (bHaveSSE)
    {
        dsp_x86_sse_add3f(fpDest, fpSrc1, fpSrc2, lCount);
    }
    else
    #endif
    {
        for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
        {
            fpDest[lLoopCntr] = fpSrc1[lLoopCntr] + fpSrc2[lLoopCntr];
        }
    }
    #endif
}


void clDSPOp::Add (double *dpDest, const double *dpSrc1, const double *dpSrc2,
    long lCount)
{
    #ifdef DSP_IPP
    ippsAdd_64f(dpSrc1, dpSrc2, dpDest, lCount);
    #else
    long lLoopCntr;

    #ifdef DSP_X86
    if (bHaveSSE)
    {
        dsp_x86_sse_add3(dpDest, dpSrc1, dpSrc2, lCount);
    }
    else
    #endif
    {
        for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
        {
            dpDest[lLoopCntr] = dpSrc1[lLoopCntr] + dpSrc2[lLoopCntr];
        }
    }
    #endif
}


void clDSPOp::Add (stpSCplx spCplxDest, const stpSCplx spCplxSrc1,
    const stpSCplx spCplxSrc2, long lCount)
{
    #ifdef DSP_IPP
    ippsAdd_32fc((Ipp32fc *) spCplxSrc1, (Ipp32fc *) spCplxSrc2,
        (Ipp32fc *) spCplxDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        CplxAdd(&spCplxDest[lLoopCntr], &spCplxSrc1[lLoopCntr],
            &spCplxSrc2[lLoopCntr]);
    }
    #endif
}


void clDSPOp::Add (stpDCplx spCplxDest, const stpDCplx spCplxSrc1,
    const stpDCplx spCplxSrc2, long lCount)
{
    #ifdef DSP_IPP
    ippsAdd_64fc((Ipp64fc *) spCplxSrc1, (Ipp64fc *) spCplxSrc2,
        (Ipp64fc *) spCplxDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        CplxAdd(&spCplxDest[lLoopCntr], &spCplxSrc1[lLoopCntr],
            &spCplxSrc2[lLoopCntr]);
    }
    #endif
}


void clDSPOp::Sub (float *fpDest, float fSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsSubC_32f_I(fSrc, fpDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        fpDest[lLoopCntr] -= fSrc;
    }
    #endif
}


void clDSPOp::Sub (double *dpDest, double dSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsSubC_64f_I(dSrc, dpDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        dpDest[lLoopCntr] -= dSrc;
    }
    #endif
}


void clDSPOp::Sub (stpSCplx spCplxDest, stSCplx sCplxSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsSubC_32fc_I(*((Ipp32fc *) &sCplxSrc), (Ipp32fc *) spCplxDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        CplxSub(&spCplxDest[lLoopCntr], &sCplxSrc);
    }
    #endif
}


void clDSPOp::Sub (stpDCplx spCplxDest, stDCplx sCplxSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsSubC_64fc_I(*((Ipp64fc *) &sCplxSrc), (Ipp64fc *) spCplxDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        CplxSub(&spCplxDest[lLoopCntr], &sCplxSrc);
    }
    #endif
}


void clDSPOp::Sub (float *fpDest, const float *fpSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsSub_32f_I(fpSrc, fpDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        fpDest[lLoopCntr] -= fpSrc[lLoopCntr];
    }
    #endif
}


void clDSPOp::Sub (double *dpDest, const double *dpSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsSub_64f_I(dpSrc, dpDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        dpDest[lLoopCntr] -= dpSrc[lLoopCntr];
    }
    #endif
}


void clDSPOp::Sub (stpSCplx spCplxDest, const stpSCplx spCplxSrc, long lCount)
{
    #ifdef DSP_IPP
    ippsSub_32fc_I((Ipp32fc *) spCplxSrc, (Ipp32fc *) spCplxDest, lCount);
    #else
    long lLoopCntr;

    for (lLoopCntr = 0L; lLoopCntr < lCount; lLoopCntr++)
    {
        CplxSub(&spCplxDest[lLoopCntr], &spCplxSrc[lLoopCntr]);
    }
    #endif
}


void clDSPOp::Sub (stpDCplx spCplxDest, const stpDCplx spCplxSrc, long lCount)
{
    #ifdef DSP_IPP