kube-gustavson-fft.hpp

Go to the documentation of this file.

#ifndef KUBE_GUSTAVSON_HPP
#define KUBE_GUSTAVSON_HPP
#include <cmath>
/**************************************************************
 * kube_fft.hpp
 *
 * Forward and inverse discrete 2D Fourier transforms.
 *
 * This software is in the public domain.
 *
 **************************************************************/

struct COMPLEX {
 float re; float im;
 // Operatoren zur einfacheren Benutzung
 COMPLEX (float real=0.0, float imag=0.0) : re(real), im(imag) {}
 COMPLEX (const COMPLEX& c) : re(c.re), im(c.im) {}
 
 COMPLEX& operator = (float val) {
   re=val;
   return *this;
 }

 COMPLEX& operator = (const COMPLEX& c) {
   re=c.re;
   im=c.im;
   return *this;
 }  
 
 COMPLEX& operator -= (const COMPLEX& c) {
   re-=c.re;
   im-=c.im;
   return *this;
 }
 
 COMPLEX operator - (const COMPLEX& c) const {
   COMPLEX erg(re-c.re, im-c.im);
   return erg;
 }
 
 COMPLEX& operator += (const COMPLEX& c) {
   re+=c.re;
   im+=c.im;
   return *this;
 }
 
 COMPLEX operator + (const COMPLEX& c) const {
   COMPLEX erg(re+c.re, im+c.im);
   return erg;
 }

 COMPLEX& operator /= (const COMPLEX& c) {
    const float r =  re * c.re + im * c.im;
    const float n =  c.re*c.re+c.im*c.im;
    im = im * c.re - re * c.im / n;
    re = r / n;
    return *this;
 }
 
 COMPLEX operator / (const COMPLEX& c) const {
   COMPLEX erg(*this);
   erg/=c;
   return erg;
 }

 bool operator != (const COMPLEX& c) const {
  
    return ((re!=c.re) || (im!=c.im)); 
 }

};

inline COMPLEX operator * (const COMPLEX& c, float mult) {
  COMPLEX erg;
  erg.re=c.re*mult;
  erg.im=c.im*mult;
  return erg;
}

inline COMPLEX operator * (float mult, const COMPLEX& c) {
  COMPLEX erg;
  erg.re=c.re*mult;
  erg.im=c.im*mult;
  return erg;
}

// Multiplikation mit einer anderen komplexen Zahl
inline COMPLEX operator *(const COMPLEX& a, const COMPLEX& b) {
  COMPLEX c;
  c.re= a.re*b.re - a.im*b.im;
  c.im= a.re*b.im + a.im*b.re;
  return c;
}  


inline COMPLEX& operator *= (COMPLEX& c, float mult) {
  c.re*=mult;
  c.im*=mult;
  return c;
}

inline float abs(const COMPLEX& c) {
  return sqrt(c.re*c.re+c.im*c.im);
}  

inline float norm(const COMPLEX& c) {
  return c.re*c.re+c.im*c.im;
}


struct DCOMPLEX {double re; double im;};

#ifndef ERROR
#define ERROR -1
#define NO_ERROR 0
#endif

#define FFT_FORWARD     0
#define FFT_INVERSE 1

#define PI      3.1415926535897932
#define TWOPI   6.2831853071795865 /* 2.0 * PI */
#define HALFPI  1.5707963267948966 /* PI / 2.0 */
#define PI8     0.392699081698724 /* PI / 8.0 */
#define RT2     1.4142135623731  /* sqrt(2.0) */
#define IRT2    0.707106781186548  /* 1.0/sqrt(2.0) */

/* Perform forward 2D transform on a COMPLEX array. */
extern int forward_fft2f(COMPLEX *array, int rows, int cols);

/* Perform inverse 2D transform on a COMPLEX array. */
extern int inverse_fft2f(COMPLEX *array, int rows, int cols);

/* Perform forward 2D transform on a DCOMPLEX array. */
extern int forward_fft2d(DCOMPLEX *array, int rows, int cols);

/* Perform inverse 2D transform on a DCOMPLEX array. */
extern int inverse_fft2d(DCOMPLEX *array, int rows, int cols);

#endif

---