FFTW_Complex.h

Go to the documentation of this file.

#ifndef FFTW_Complex_H
#define FFTW_Complex_H



#include "FFTW_Object.h"

#include "OMP_Thread.h"


DEFINE_SPTR(FFTW_Complex);
class FFTW_Complex : public virtual FFTW_Object {
    SP_OBJECT(FFTW_Complex);
    // ATTRIBUTES


private:
    //complex value
    tFFTWComplex mValue;

    //real value
    mutable tFFTWReal mWork;
    // ASSOCIATIONS
    
    
public:
    // METHODS
    
    // CONSTRUCTORS 
    
    FFTW_Complex(void) {
        mValue[0]=0;
        mValue[0]=0;
    }
    FFTW_Complex(const tFFTWReal& a,const tFFTWReal& b) {
        mValue[0]=a;
        mValue[1]=b;
    }
    FFTW_Complex(const tFFTWComplex& f) {
        mValue[0]=f[0];
        mValue[1]=f[1];
    }
    FFTW_Complex(const FFTW_Complex& f) {
        mValue[0]=f[0];
        mValue[1]=f[1];
    }
   
    // DESTRUCTORS  
    virtual ~FFTW_Complex(void) {
    }
    
    // -----------------------
    // shared pointer operator
    // ------------------------
public:
  
    // ----------------
    // New constructors
    // ----------------
public:
    static inline SP::FFTW_Complex New() {
        SP::FFTW_Complex p(new FFTW_Complex(),
                           CORE_Object::Delete());
        p->setThis(p);
        return p;
    };
    
public:
    // operators
    // =========
    inline tFFTWReal& operator[](const tUSInt& i)  {
        return mValue[i];
    };
    inline const tFFTWReal& operator[](const tUSInt& i)  const {
        return mValue[i];
    };
  
    inline FFTW_Complex & operator=(const FFTW_Complex&  s)  {
        mValue[0]=s[0];
        mValue[1]=s[1];
        return (*this);
    };
    inline FFTW_Complex& operator=(const tFFTWComplex&  s)  {
        mValue[0]=s[0];
        mValue[1]=s[1];
        return (*this);
    };
   
    inline FFTW_Complex& operator=(const tFFTWReal&  s)  {
        mValue[0]=s;
        mValue[1]=0;
        return (*this);
    };
   
    inline FFTW_Complex& operator*=(const tFFTWReal&  s)  {
        mValue[0]*=s;
        mValue[1]*=s;
        return (*this);
    };
    inline FFTW_Complex& operator*=(const FFTW_Complex&  s)  {
        mWork=mValue[0]*s[0]-mValue[1]*s[1];
        mValue[1]=mValue[0]*s[1]+mValue[1]*s[0];
        mValue[0]=mWork;
        return (*this);
    };
    
    inline FFTW_Complex& operator*=(const tFFTWComplex&  s)  {
        mWork=mValue[0]*s[0]-mValue[1]*s[1];
        mValue[1]=mValue[0]*s[1]+mValue[1]*s[0];
        mValue[0]=mWork;
        return (*this);
    };
    
    inline FFTW_Complex& operator/=(const tFFTWReal&  s)  {
        if (s==0) throw CORE_Exception("math/signam","FFTW_Complex/=","division by 0");
        mValue[0]/=s;
        mValue[1]/=s;
        return (*this);
    };
    inline FFTW_Complex& operator/=(const FFTW_Complex&  s)  {
        mWork=mValue[0]*s[0]+mValue[1]*s[1];
        mValue[1]=(-mValue[0]*s[1]+mValue[1]*s[0]);
        mValue[0]=mWork;
        mWork=s.module2();
        if (mWork==0) throw CORE_Exception("math/signam","FFTW_Complex/=","division by 0");
        mValue[0]/=mWork;
        mValue[1]/=mWork;
        return (*this);
    };
    inline FFTW_Complex& operator+=(const FFTW_Complex&  s)  {
        mValue[0]+=s[0];
        mValue[1]+=s[1];
        return (*this);
    };
    inline FFTW_Complex& operator+=(const tFFTWComplex&  s)  {
        mValue[0]+=s[0];
        mValue[1]+=s[1];
        return (*this);
    };
    inline FFTW_Complex& operator-=(const FFTW_Complex&  s)  {
        mValue[0]-=s[0];
        mValue[1]-=s[1];
        return (*this);
    };
    inline FFTW_Complex& operator-=(const tFFTWComplex&  s)  {
        mValue[0]-=s[0];
        mValue[1]-=s[1];
        return (*this);
    };
    // SET methods

    inline void setValue(const tFFTWReal& x,const tFFTWReal& y) {
        mValue[0]=x;
        mValue[1]=y;
    }
    
    // GET methods

    inline const tFFTWComplex& getValue() const {
        return mValue;
    }
    inline tFFTWComplex& getValue()  {
        return mValue;
    }

    inline void conjugate() {
        mValue[1]*=-1;
    }
    inline tFFTWReal module2() const {
        return module2(mValue);
    }
    inline tFFTWReal module() const {
        return module(mValue);
    }
    inline static tFFTWReal module2(const tFFTWComplex& f)  {
        return f[0]*f[0]+f[1]*f[1];
    }
    inline static tFFTWReal module(const tFFTWComplex& f) {
        return sqrt(module2(f));
    }
    
    /* \brief  complex product a*=b;
     * @param a: the a value
     * @param b: the b value  
     */
    inline static void product(tFFTWComplex&  a,const tFFTWComplex& b) {
        tFFTWReal temp;
        product(a,b,temp);
    }
    /* \brief  complex product a*=b;
     * @param a: the a value
     * @param b: the b value
     */
    inline static void product(tFFTWComplex&  a,const tFFTWComplex& b,tFFTWReal& wrk) {
        product(a,b,a,wrk);
    }
    
    /* \brief  complex product r=a*b;
     * @param a: the a value
     * @param b: the b value
     * @param r: the result value
     * @param temp: work value
     */
    inline static void product(const tFFTWComplex&  a,const tFFTWComplex& b,tFFTWComplex& r,tFFTWReal& wrk) {
        wrk=a[0]*b[0]-a[1]*b[1];
        r[1]=a[0]*b[1]+a[1]*b[0];
        r[0]=wrk;
    }
    /* \brief  complex product r=a*b
     * @param a: the a value
     * @param b: the b value
     * @param r: the result value
     */
    inline static void product(const tFFTWComplex&  a,const tFFTWComplex& b,tFFTWComplex& r) {
        ASSERT_IN(&b!=&r);
        ASSERT_IN(&a!=&r);
        r[0]=a[0]*b[0]-a[1]*b[1];
        r[1]=a[0]*b[1]+a[1]*b[0];
    }
    /* \brief  complex product r=a*b;
     * @param a: the a value
     * @param b: the b value
     * @param r: the result value
     */
    inline static void product(const FFTW_Complex&  a,const FFTW_Complex& b,FFTW_Complex& r) {
        product(a.getValue(),b.getValue(),r.getValue());
    }
    /* \brief  complex product r=(a+ib)*(x+i*y)
     * @param a: real part
     * @param b: imaginary part
     * @param x: real part
     * @param y: imaginary part
     * @param r: result 
     */
    inline static void product(const tFFTWReal&  a,const tFFTWReal& b,const tFFTWReal& x,const tFFTWReal& y,tFFTWComplex& r) {
        r[0]=a*x-b*y;
        r[1]=a*y+b*x;
    }
    /* \brief  complex add a+=b
     * @param a: the a value
     * @param b: the b value
     */
    inline static void add(tFFTWComplex&  a,const tFFTWComplex& b) {
        a[0]+=b[0];
        a[1]+=b[1];
    }
    /* \brief  complex add  r=a+b;
     * @param a: the a value
     * @param b: the b value
     * @param r: the result value
     */
    inline static void add(const tFFTWComplex&  a,const tFFTWComplex& b,tFFTWComplex& r) {
        r[0]=a[0]+b[0];
        r[1]=a[0]+b[1];
    }
    /* \brief  complex sub a-=b;
     * @param a: the a value
     * @param b: the b value
     */
    inline static void sub(tFFTWComplex&  a,const tFFTWComplex& b) {
        a[0]-=b[0];
        a[1]-=b[1];
    }
    /* \brief  complex sub : r=a-b;
     * @param a: the a value
     * @param b : the b value
     * @param r: the result value
     */
    inline static void sub(const tFFTWComplex&  a,const tFFTWComplex& b,tFFTWComplex& r) {
        r[0]=a[0]-b[0];
        r[1]=a[1]-b[1];
    }
    /* \brief  swap a & b
     * @param a: the value to swap
     * @param b: the value to swap
     */
    inline static void swap(tFFTWComplex& a,tFFTWComplex& b) {
        std::swap(a[0],b[0]);
        std::swap(a[1],b[1]);
    }

    inline static tFFTWReal norm2(const tFFTWComplex&  a,const tFFTWComplex& b) {
       
        return (a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1]);
    }
    inline static tFFTWReal norm(const tFFTWComplex&  a,const tFFTWComplex& b) {
        
        return sqrt(norm2(a,b));
    }



    virtual tString toString() const {
        return toString(mValue);
    }
    inline static tString toString(const tFFTWComplex& f) {
        tString ret="";
        ret+=CORE_Real::toString(f[0]);
        if (f[1]>0) {
            ret+="+"+CORE_Real::toString(f[1])+"i";
        } else if (f[1]<0) {
            ret+=CORE_Real::toString(f[1])+"i";
        }
        return ret;
    }
    inline static tString toString(const fftw_complex& f,const tUSInt& nDigits) {
        tString ret="";
        ret+=CORE_Real::toString(f[0],nDigits);
        if (f[1]>0) {
            ret+="+"+CORE_Real::toString(f[1],nDigits)+"i";
        } else if (f[1]<0) {
            ret+=CORE_Real::toString(f[1],nDigits)+"i";
        }
        return ret;
    }
    
    
    
    
    
};


#endif