Implement 2-dimensional Fast Fourier Transformation (#1731)

Jürgen Knödlseder
1 parent 79d739e6
Showing 8 changed files with 241 additions and 72 deletions
include/GFft.hpp
src/model/GModelSpectralExpInvPlaw.cpp
src/numerics/GFft.cpp
src/numerics/GFftWavetable.cpp
src/numerics/GNdarray.cpp
src/vo/GVOTable.cpp
test/test_GNumerics.cpp
test/test_GNumerics.hpp
@@ -66,7 +66,9 @@ public:
  
     // FFT access operators
     std::complex<double>&       operator()(const int& ix);
+    std::complex<double>&       operator()(const int& ix, const int& iy);
     const std::complex<double>& operator()(const int& ix) const;
+    const std::complex<double>& operator()(const int& ix, const int& iy) const;
  
     // Operators
     GFft& operator=(const GFft& fft);
@@ -197,7 +199,7 @@ std::string GFft::classname(void) const
 /***********************************************************************//**
  * @brief 1-dimensional FFT element access operator
  *
- * @param[in] ix Element index [0,...,size()-1].
+ * @param[in] ix Element index [0,...,shape(0)-1].
  * @return Reference to FFT element.
  ***************************************************************************/
 inline
@@ -209,9 +211,24 @@ std::complex&lt;double&gt;&amp; GFft::operator()(const int&amp; ix)
  
  
 /***********************************************************************//**
+ * @brief 2-dimensional FFT element access operator
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @return Reference to FFT element.
+ ***************************************************************************/
+inline
+std::complex<double>& GFft::operator()(const int& ix, const int& iy)
+{
+	// Return array element
+    return *(m_data + ix + m_strides[1]*iy);
+}
+
+
+/***********************************************************************//**
  * @brief 1-dimensional FFT element access operator (const variant)
  *
- * @param[in] ix Element index [0,...,size()-1].
+ * @param[in] ix Element index [0,...,shape(0)-1].
  * @return Const reference to FFT element.
  ***************************************************************************/
 inline
@@ -223,6 +240,23 @@ const std::complex&lt;double&gt;&amp; GFft::operator()(const int&amp; ix) const
  
  
 /***********************************************************************//**
+ * @brief 2-dimensional FFT element access operator (const variant)
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @return Const reference to FFT element.
+ ***************************************************************************/
+inline
+const std::complex<double>& GFft::operator()(const int& ix, const int& iy) const
+{
+	// Return array element
+    return *(m_data + ix + m_strides[1]*iy);
+}
+
+
+
+
+/***********************************************************************//**
  * @brief Return dimension of Fast Fourier Transformation
  *
  * @return Dimension of Fast Fourier Transformation.
@@ -85,7 +85,7 @@ GModelSpectralExpInvPlaw::GModelSpectralExpInvPlaw(void) : GModelSpectral()
  * @param[in] prefactor Name of prefactor parameter.
  * @param[in] index Name of index parameter.
  * @param[in] pivot Name of pivot parameter.
- * @param[in] cutoff Name of cutoff parameter.
+ * @param[in] lambda Name of cutoff parameter.
  ***************************************************************************/
 GModelSpectralExpInvPlaw::GModelSpectralExpInvPlaw(const std::string& type,
                                                    const std::string& prefactor,
@@ -35,6 +35,8 @@
 #include "GNdarray.hpp"
  
 /* __ Method name definitions ____________________________________________ */
+#define G_FORWARD                                  "GFft::forward(GNdarray&)"
+#define G_BACKWARD                                         "GFft::backward()"
  
  
 /*==========================================================================
@@ -277,7 +279,27 @@ GFft* GFft::clone(void) const
  *
  * @param[in] array N-dimensional array.
  *
- * @todo Method should support N-dim arrays, only supports 1-dim so far
+ * GException::feature_not_implemented
+ *             Array dimension is not supported.
+ *
+ * Performs discret Fast Fourier forward transformation for 1- and
+ * 2-dimensional arrays.
+ *
+ * For one dimensional arrays of length \f$N\f$ the transformation is given
+ * by
+ *
+ * \f[
+ *    X_k = \sum_{n=0}^{N-1} z_n e^{-i 2 \pi \frac{k n}{N}}
+ * \f]
+ *
+ * and for two dimensional arrays of length \f$ M \times N\f$ the
+ * transformation is given by
+ *
+ * \f[
+ *    X_{k,l} = \sum_{n=0}^{N-1} \left(
+ *              \sum_{m=0}^{M-1} z_{m,n} e^{-i 2 \pi \frac{m k}{M}} \right)
+ *              e^{-i 2 \pi \frac{n l}{N}}
+ * \f]
  ***************************************************************************/
 void GFft::forward(const GNdarray& array)
 {
@@ -286,26 +308,36 @@ void GFft::forward(const GNdarray&amp; array)
  
     // Continue only if there are elements in the data member
     if (m_size > 0) {
-/*
-        // Loop over all dimensions of array
-        for (int d = 0; d < dim(); ++d) {
  
-            // Loop over all dimensions of array
-            for (int i = 0; i < dim(); ++i) {
+        // Case A: array is 1-dimensional
+        if (array.dim() == 1) {
+            transform(m_data, 1, m_shape[0], m_wavetable[0], true);
+        }
  
-                // Skip
-                if (d == i) {
-                    continue;
-                }
+        // Case B: array is 2-dimensional
+        else if (array.dim() == 2) {
  
-                // Loop over all
-                int n      = shape()[d];
-                int stride = strides()[d];
-                int istart =
-*/
+            // Loop over rows
+            for (int row = 0; row < m_shape[1]; ++row) {
+                transform(m_data + row*m_strides[1], m_strides[0], m_shape[0],
+                          m_wavetable[0], true);
+            }
  
-        // Perform FFT (is only correct for 1D for now)
-        transform(m_data, 1, m_shape[0], m_wavetable[0], true);
+            // Loop over columns
+            for (int col = 0; col < m_shape[0]; ++col) {
+                transform(m_data + col*m_strides[0], m_strides[1], m_shape[1],
+                          m_wavetable[1], true);
+            }
+
+        }
+
+        // ... otherwise throw an exception
+        else {
+            std::string msg = "Fast Fourier Transformation for "+
+                              gammalib::str(array.dim())+"-dimensional array "
+                              "not implemented.";
+            throw GException::feature_not_implemented(G_FORWARD, msg);
+        }
  
     } // endif: there were elements in the FFT data member
  
@@ -319,7 +351,28 @@ void GFft::forward(const GNdarray&amp; array)
  *
  * @return N-dimensional array.
  *
- * @todo Method should support N-dim arrays, only supports 1-dim so far
+ * GException::feature_not_implemented
+ *             Array dimension is not supported.
+ *
+ * Performs the inverse (or backward) discret Fast Fourier forward
+ * transformation for 1- and 2-dimensional arrays.
+ *
+ * For one dimensional arrays of length \f$N\f$ the transformation is given
+ * by
+ *
+ * \f[
+ *    z_k = \frac{1}{N} \sum_{n=0}^{N-1} X_n e^{i 2 \pi \frac{k n}{N}}
+ * \f]
+ *
+ * and for two dimensional arrays of length \f$ M \times N\f$ the
+ * transformation is given by
+ *
+ * \f[
+ *    z_{k,l} = \frac{1}{N} \sum_{n=0}^{N-1} \left(
+ *              \frac{1}{M} \sum_{m=0}^{M-1} X_{m,n}
+ *                          e^{i 2 \pi \frac{m k}{M}} \right)
+ *                          e^{i 2 \pi \frac{n l}{N}}
+ * \f]
  ***************************************************************************/
 GNdarray GFft::backward(void) const
 {
@@ -327,24 +380,66 @@ GNdarray GFft::backward(void) const
     GNdarray array(shape());
  
     // Continue only if there are elements in the array
-    if (array.size() > 0) {
+    if (m_size > 0) {
  
         // Create copy of complex array
-        std::complex<double>* data = new std::complex<double>[array.size()];
-        for (int i = 0; i < array.size(); ++i) {
+        std::complex<double>* data = new std::complex<double>[m_size];
+        for (int i = 0; i < m_size; ++i) {
             *(data+i) = *(m_data+i);
         }
  
-        // Perform FFT (circumvent const correctness)
-        const_cast<GFft*>(this)->transform(data, 1, m_shape[0], m_wavetable[0],
-                                           false);
+        // Case A: array is 1-dimensional
+        if (array.dim() == 1) {
+
+            // Perform backward transformation (circumvent const correctness)
+            const_cast<GFft*>(this)->transform(data, 1, m_shape[0], m_wavetable[0], false);
+
+            // Get normalisation factor
+            const double norm = 1.0 / double(m_shape[0]);
+
+            // Extract 1-dimensional array and normalise with 1/n
+            for (int i = 0; i < m_size; ++i) {
+                array(i) = (data+i)->real() * norm;
+            }
+
+        }
+
+        // Case B: array is 2-dimensional
+        else if (array.dim() == 2) {
  
-        // Get normalisation factor
-        const double norm  = 1.0 / double(array.size());
+            // Loop over rows
+            for (int row = 0; row < m_shape[1]; ++row) {
+                const_cast<GFft*>(this)->transform(data + row*m_strides[1], m_strides[0], m_shape[0],
+                          m_wavetable[0], false);
+            }
+
+            // Loop over columns
+            for (int col = 0; col < m_shape[0]; ++col) {
+                const_cast<GFft*>(this)->transform(data + col*m_strides[0], m_strides[1], m_shape[1],
+                          m_wavetable[1], false);
+            }
+
+            // Get normalisation factor
+            const double norm = 1.0 / double(m_shape[0]*m_shape[1]);
+
+            // Extract 2-dimensional array and normalise with 1/(n*m)
+            for (int i = 0; i < m_size; ++i) {
+                array(i) = (data+i)->real() * norm;
+            }
  
-        // Extract N-dimensional array and normalise with 1/n
-        for (int i = 0; i < array.size(); ++i) {
-            array(i) = (data+i)->real() * norm;
+        }
+
+        // ... otherwise throw an exception
+        else {
+
+            // Delete copy of complex array
+            delete [] data;
+
+            // Throw exception
+            std::string msg = "Fast Fourier Transformation for "+
+                              gammalib::str(array.dim())+"-dimensional array "
+                              "not implemented.";
+            throw GException::feature_not_implemented(G_BACKWARD, msg);
         }
  
         // Delete copy of complex array
@@ -501,15 +596,9 @@ void GFft::set_data(const GNdarray&amp; array)
             *(m_data+i) = std::complex<double>(array(i), 0.0);
         }
  
-        // Save shape
-        m_shape = array.shape();
-
-        // Set strides
-        int stride = 1;
-        for (int i = 0; i < m_shape.size(); ++i) {
-            m_strides.push_back(stride);
-            stride *= m_shape[i];
-        }
+        // Save shape and strides
+        m_shape   = array.shape();
+        m_strides = array.strides();
  
         // Set trigonometric coefficients
         for (int i = 0; i < m_shape.size(); ++i) {
@@ -1319,30 +1408,18 @@ void GFft::factorn(std::complex&lt;double&gt;* in,
             for (int e1 = 1; e1 < factor; ++e1) {
  
                 // Get x
-                double x_real = (in+istride*(i + e1 * m))->real();
-                double x_imag = (in+istride*(i + e1 * m))->imag();
+                std::complex<double> x = *(in+istride*(i + e1 * m));
  
-                // Get w
-                double w_real, w_imag;
-             
-                // Compute indices
+                // Compute index
                 int twiddle = index + (e1-1)*q + k-1;
  
-                // Set trigonometric coefficients for forward transform
-                if (sign == -1) {
-                    w_real = wavetable[twiddle].real();
-                    w_imag = wavetable[twiddle].imag();
-                }
-                // ... otherwise set trigonometric coefficients for backward
-                // tranform: w -> conjugate(w)
-                else {
-                    w_real =  wavetable[twiddle].real();
-                    w_imag = -wavetable[twiddle].imag();
-                }
+                // Get w
+                std::complex<double> w = (sign == -1)
+                                         ? wavetable[twiddle]
+                                         : std::conj(wavetable[twiddle]);
  
-                //
-                *(out+ostride*(j + e1 * p_1)) = std::complex<double>(w_real * x_real - w_imag * x_imag,
-                                                                     w_real * x_imag + w_imag * x_real);
+                // Compute out = w * x
+                *(out+ostride*(j + e1 * p_1)) = w * x;
  
             } // endfor: e1
  
@@ -383,7 +383,7 @@ void GFftWavetable::set_members(const int&amp; n)
  *
  * @exception GException::invalid_argument
  *            Non positive array length specified
- * @expection GException::invalid_value
+ * @exception GException::invalid_value
  *            Product of factorisation factors are not equal to array length
  *
  * Computes the factorisation
@@ -777,6 +777,7 @@ bool GNdarray::has_same_shape(const GNdarray&amp; array) const
 /***********************************************************************//**
  * @brief Throw exception if array shapes differ
  *
+ * @param[in] method Method that throws exception.
  * @param[in] array Array.
  ***************************************************************************/
 void GNdarray::require_same_shape(const std::string& method,
@@ -120,7 +120,7 @@ GVOTable::~GVOTable(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] VO table.
+ * @param[in] votable VO table.
  * @return VO table.
  ***************************************************************************/
 GVOTable& GVOTable::operator=(const GVOTable& votable)
@@ -282,7 +282,7 @@ void TestGNumerics::test_ndarray(void)
  ***************************************************************************/
 void TestGNumerics::test_fft(void)
 {
-    // Loop over array lengths
+    // Test 1-dimensional FFT
     for (int n = 2; n < 12; ++n) {
  
         // Skip n=8,9,10 since they can be factorised
@@ -303,25 +303,32 @@ void TestGNumerics::test_fft(void)
         GNdarray back = fft.backward();
  
         // Check FFT results
-        check_fft(array, fft, back);
+        check_fft1(array, fft, back);
  
     } // endfor: looped over array lengths
  
+    // Test 2-dimensional FFT
+    GNdarray array2(8,9);
+    array2(4,4) = 1.0;
+    GFft     fft2(array2);
+    GNdarray back2 = fft2.backward();
+    check_fft2(array2, fft2, back2);
+
     // Exit test
     return;
 }
  
  
 /***********************************************************************//**
- * @brief Check Fast Fourier Transform
+ * @brief Check 1-dimensional Fast Fourier Transform
  *
  * @param[in] array Input array.
  * @param[in] fft Fast Fourier Transform of input array.
  * @param[in] back Backward transform of FFT (should be equal to input array)
  ***************************************************************************/
-void TestGNumerics::check_fft(const GNdarray& array,
-                              const GFft&     fft,
-                              const GNdarray& back)
+void TestGNumerics::check_fft1(const GNdarray& array,
+                               const GFft&     fft,
+                               const GNdarray& back)
 {
     // Get normalisation factor
     double norm = 1.0 / double(array.size());
@@ -350,6 +357,53 @@ void TestGNumerics::check_fft(const GNdarray&amp; array,
  
  
 /***********************************************************************//**
+ * @brief Check 2-dimensional Fast Fourier Transform
+ *
+ * @param[in] array Input array.
+ * @param[in] fft Fast Fourier Transform of input array.
+ * @param[in] back Backward transform of FFT (should be equal to input array)
+ ***************************************************************************/
+void TestGNumerics::check_fft2(const GNdarray& array,
+                               const GFft&     fft,
+                               const GNdarray& back)
+{
+    // Get normalisation factor
+    double normM = 1.0 / double(array.shape()[0]);
+    double normN = 1.0 / double(array.shape()[1]);
+
+    // Compute reference value and compare with FFT result
+    for (int k = 0; k < array.shape()[0]; ++k) {
+        for (int l = 0; l < array.shape()[1]; ++l) {
+            double ref_real = 0.0;
+            double ref_imag = 0.0;
+            for (int n = 0; n < array.shape()[1]; ++n) {
+                double sum_real = 0.0;
+                double sum_imag = 0.0;
+                for (int m = 0; m < array.shape()[0]; ++m) {
+                    double x = -gammalib::twopi * m * k * normM;
+                    sum_real += array(m,n) * std::cos(x);
+                    sum_imag += array(m,n) * std::sin(x);
+                }
+                double y = -gammalib::twopi * n * l * normN;
+                ref_real += (sum_real * std::cos(y) - sum_imag * std::sin(y));
+                ref_imag += (sum_real * std::sin(y) + sum_imag * std::cos(y));
+            }
+            std::complex<double> ref(ref_real, ref_imag);
+            test_value(fft(k,l), ref);
+        }
+    }
+
+    // Test backward transformation
+    for (int k = 0; k < array.size(); ++k) {
+        test_value(back(k), array(k));
+    }
+
+    // Exit test
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Test model parameter handling
  ***************************************************************************/
 void TestGNumerics::test_integral(void)
@@ -70,9 +70,12 @@ public:
     virtual std::string    classname(void) const { return "TestGNumerics"; }
     void                   test_ndarray(void);
     void                   test_fft(void);
-    void                   check_fft(const GNdarray& array,
-                                     const GFft&     fft,
-                                     const GNdarray& back);
+    void                   check_fft1(const GNdarray& array,
+                                      const GFft&     fft,
+                                      const GNdarray& back);
+    void                   check_fft2(const GNdarray& array,
+                                      const GFft&     fft,
+                                      const GNdarray& back);
     void                   test_integral(void);
     void                   test_romberg_integration(void);
     void                   test_adaptive_simpson_integration(void);