|| || || || || || || ||

---

MakeGaborIIR1d.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_MakeGaborIIR1d_h
#define Impala_Core_Array_MakeGaborIIR1d_h

#include "Util/Math.h"
#include "Core/Array/Arrays.h"

namespace Impala
{
namespace Core
{
namespace Array
{


inline Array2dComplex64*
MakeGaborIIR1d(double sigma, double omega0,
               Complex64 &borderIn, Complex64 &borderOut)
{
    int size = 7;
    double _re[7];
    double _im[7];

    // do the 'recipe' in the article

    // initial values
    double m0 = 1.16680, m1 = 1.10783, m2 = 1.40586 ;
    double m1sq = m1*m1, m2sq = m2*m2 ;

    // calculate q
    double q ;
    if(sigma < 3.556)
        q = -0.2568 + 0.5784 * sigma + 0.0561 * sigma * sigma ;
    else
        q = 2.5091 + 0.9804 * (sigma - 3.556) ;

    double qsq = q*q ;

    // calculate scale, and b[0,1,2,3]
    double scale = (m0 + q) * (m1sq + m2sq + 2*m1*q + qsq) ;
    double b0 = 1 ;
    double b1 = -q * (2*m0*m1 + m1sq + m2sq + (2*m0 + 4*m1) * q + 3*qsq) / scale ;
    double b2 = qsq * (m0 + 2*m1 + 3*q) / scale ;
    double b3 = - qsq * q / scale ;

    // calculate B, without the ^2, this because horus uses the middle value of
    // the filter both for left and right parts of the recursion
    double B = (m0 * (m1sq + m2sq))/scale  ;

    // create filter, 
    // yet, as we use addAssign, and the article uses subAssign,
    // we use the negative values 
    _re[0] = -b3 * cos(3*omega0)    ;  _im[0] =  -b3 * sin(3*omega0) ;
    _re[1] = -b2 * cos(2*omega0)    ;  _im[1] =  -b2 * sin(2*omega0) ;
    _re[2] = -b1 * cos(omega0)      ;  _im[2] =  -b1 * sin(omega0) ;
    _re[3] = B                      ;  _im[3] = 0 ;
    _re[4] = _re[2]                 ;  _im[4] =  -_im[2] ;
    _re[5] = _re[1]                 ;  _im[5] =  -_im[1] ;
    _re[6] = _re[0]                 ;  _im[6] =  -_im[0] ;


    // create border coefficients as done in the article,
    double leftReal =       1.0 - _re[2] - _re[1] - _re[0] ;
    double leftImaginary =      - _im[2] - _im[1] - _im[0] ;

    double rightReal =      1.0 - _re[4] - _re[5] - _re[6] ;
    double rightImaginary =     - _im[4] - _im[5] - _im[6] ;

    borderIn = Complex64(leftReal,leftImaginary) ;
    borderOut = Complex64(rightReal,rightImaginary) ;

    Complex64 complexB(B,0);
    Complex64 one(1.0,0.0) ;  

    // normalize the 'in' filter also for B, as horus uses B also for the 'in' part.
    borderIn = complexB / borderIn ;
    borderOut = complexB / borderOut  ;

    //Complex* result = new Complex[7];
    Array2dComplex64* result = ArrayCreate<Array2dComplex64>(7, 1);
    Real64* ptr = result->CPB(0, 0);
    for (int i=0 ; i<7 ; i++)
    {
        //result[i] = Complex(_re[i], _im[i]);
        *ptr++ = _re[i];
        *ptr++ = _im[i];
    }
    return result;
}

} // namespace Array
} // namespace Core
} // namespace Impala

#endif

---