Changeset 314

Timestamp:

Mar 29, 2004, 7:40:01 PM (22 years ago)

Author:

Paul Price

Message:

Major hack of the FFT APIs, incorporating Nick's comments: Defined a
FFT type, on which all the FFT operations act. This insulates us from
the details of the implementation (e.g., FFTW requires a "plan";
different FFT libraries store the transformed data differently).
Added the functionality Nick suggested (filtering, cross-correlation).

Location:

trunk

Files:

• archive/pslib/include/psFFT.h (modified) (2 diffs)
• doc/pslib/psLibSDRS.tex (modified) (2 diffs)

trunk/archive/pslib/include/psFFT.h

@@ -1 +1 @@
 #if !defined(PS_FFT_H)
 #define PS_FFT_H
+
+#include <fftw3.h>
 
 /** \file psFFT.h
@@ -12 +14 @@
  */
 
-/** Return Fourier Transform of an array */
-psComplexArray *
-psRealFFT(psComplexArray *restrict out, //!< Output array to be returned; may be NULL
-          const psFloatArray *restrict myArray //!< Input array
+/** Details on FFT implementation (private).  Example shown is for FFTW */
+typedef struct {
+    fftw_plan plan;                     //!< FFTW plan on how to do the FFT
+    char *filename;                     //!< File name for FFTW plan
+} p_psFFTDetails;
+
+/** Fast Fourier Transform */
+typedef struct {
+    p_psFFTDetails *details;            //!< Details on FFT implementation (private)
+    int nx, ny;                         //!< Size in x and y
+    float **real;                       //!< Data in real space
+    void *fourier;                      //!< Data in fourier space; implementation dependent
+} psFFT;
+
+/** Constructor */
+psFFT *
+psFFTAlloc(psImage *image               //!< Image to transform
+           );
+
+/** Constructor for 1D case */
+psFFT *
+psFFTAlloc1D(const psFloatArray *arr    //!< Array to transform
+             );
+
+/** Destructor. Returns the data in the real space as an image. */
+psImage *
+psFFTFree(psImage *out,                 //!< Image to write the data to, or NULL
+          psFFT *restrict fft           //!< FFT to destroy
           );
 
-/** Return [inverse?] Fourier Transform of a complex array.  May NOT be done in-place (restrict-ed) */
-psComplexArray *
-psComplexFFT(psComplexArray *restrict out, //!< Output array to be returned; may be NULL
-             const psComplexArray *restrict myArray, //!< Input array
-             int sign                   //!< +1 or -1 to indicate direction of FT
-             );
+/** Forward FFT: from real to fourier space */
+psFFT *
+psFFTForwardTransform(psFFT *fft        //!< FFT to apply (input and output)
+                      );
 
-/** Return Power spectrum of a array */
-psFloatArray *
-psPowerSpec(psFloatArray *restrict out, //!< Output array to be returned
-            const psFloatArray *restrict myArray //!< Input array
+/** Reverse FFT: from fourier to real space */
+psFFT *
+psFFTReverseTransform(psFFT *fft        //!< FFT to apply (input and output)
+                      );
+
+/** Apply filter function in fourier space */
+psFFT *
+psFFTFilter(psFFT *fft,                 //!< FFT to use (input and output)
+            float (*filterFunc)(int kx, int ky) //!< External filter function
             );
+
+/** Apply complex filter function */
+psFFT *
+psFFTFilterComplex(psFFT *fft,          //!< FFT to use (input and output)
+                   float (*realFilterFunc)(int kx, int ky), //!< External filter function, real part
+                   float (*imagFilterFunc)(int kx, int ky) //!< External filter function, imaginary part
+                   );
+
+/** Calculate cross-correlation function */
+psFFT *
+psFFTCrossCorrelate(psFFT *out          //!< Output FFT (or NULL)
+                    const psFFT *fft1, const psFFT *fft2 //!< FFTs to use in cross-correlation
+                    );
+
+/** Calculate power spectrum */
+psFFT *
+psFFTPowerSpec(psFFT *fft               //!< FFT to use (input and output)
+               );
+
+/** Multiply two Fourier transforms, as for convolution */
+psFFT *
+psFFTMultiplyFT(psFFT *out,             //!< Output FFT (or NULL)
+                const psFFT *fft1, const psFFT *fft2 //!< FFTs to multiply
+                );
+
+/* Convert the real data in the FFT struct to an image again */
+psImage *
+psFFTGetImage(psImage *out,             //!< Image to write to (or NULL)
+              const psFFT *fft          //!< FFT to get image from
+              );
+
+/** Convert the Fourier transform data in the FFT struct to an image of complex numbers */
+psImage *
+psFFTGetFT(psImage *out,                //!< Image to write to (or NULL)
+           const psFFT *fft             //!< FFT to get Fourier transform from
+           );
 
 /* \} */ // End of MathGroup Functions
 
 #endif
+
+

trunk/doc/pslib/psLibSDRS.tex

@@ -1 @@
-%%% $Id: psLibSDRS.tex,v 1.13 2004-03-29 22:13:13 price Exp $
+%%% $Id: psLibSDRS.tex,v 1.14 2004-03-30 05:40:01 price Exp $
 \documentclass[panstarrs]{panstarrs}
 %\documentclass[panstarrs]{panstarrs}
@@ -1437 +1437 @@
 
 We require the ability to calculate the (fast) fourier transforms of
-both floating-point and complex floating-point type arrays, as well as
-the power spectrum.  The APIs are:
-
-\begin{verbatim}
-/** Return Fourier Transform of an array */
-psComplexArray *
-psRealFFT(psComplexArray *restrict out, //!< Output array to be returned; may be NULL
-          const psFloatArray *restrict myArray //!< Input array
-          );
-\end{verbatim}
-
-\begin{verbatim}
-/** Return [inverse?] Fourier Transform of a complex array.  May NOT be done in-place (restrict-ed) */
-psComplexArray *
-psComplexFFT(psComplexArray *restrict out, //!< Output array to be returned; may be NULL
-             const psComplexArray *restrict myArray, //!< Input array
-             int sign                   //!< +1 or -1 to indicate direction of FT
-             );
-\end{verbatim}
-
-\begin{verbatim}
-/** Return Power spectrum of a array */
-psFloatArray *
-psPowerSpec(psFloatArray *restrict out, //!< Output array to be returned
-            const psFloatArray *restrict myArray //!< Input array
-            );
+floating-point two-dimensional arrays (including one-dimensional
+arrays as a limiting case), and use these to perform filtering, and
+calculate the power spectrum and cross-correlation.  We expect that
+these will be implemented through wrapping an external library.  We
+define a FFT type:
+
+\begin{verbatim}
+/** Fast Fourier Transform */
+typedef struct {
+    p_psFFTDetails *details;            //!< Details on FFT implementation (private)
+    int nx, ny;                         //!< Size in x and y
+    float **real;                       //!< Data in real space
+    void *fourier;                      //!< Data in fourier space; implementation dependent
+} psFFT;
+\end{verbatim}
+
+The \code{p_psFFTDetails *details} entry allows us to wrap the
+external library in the event that it carries information around from
+one transform to the next (as is the case for FFTW).  The intent is
+that every function that employs a FFT will act on a \code{psFFT}.  A
+user will create a \code{psFFT} by stuffing it with the real data,
+perform operations in Fourier space using the \code{psFFT}, and then
+extract the output as an image when done.
+
+We specify two constructors --- one for use with images, and the other
+for use with a one-dimensional array.  The destructor returns an image
+that includes the real space data.
+
+\begin{verbatim}
+/** Constructor */
+psFFT *
+psFFTAlloc(psImage *image               //!< Image to transform
+           );
+
+/** Constructor for 1D case */
+psFFT *
+psFFTAlloc1D(const psFloatArray *arr    //!< Array to transform
+             );
+
+/** Destructor. Returns the data in the real space as an image. */
+psImage *
+psFFTFree(psImage *out,                 //!< Image to write the data to, or NULL
+          psFFT *restrict fft           //!< FFT to destroy
+          );
+\end{verbatim}
+
+The forward and reverse Fourier transforms may be calculated:
+
+\begin{verbatim}
+/** Forward FFT: from real to fourier space */
+psFFT *
+psFFTForwardTransform(psFFT *fft        //!< FFT to apply
+                      );
+
+/** Reverse FFT: from fourier to real space */
+psFFT *
+psFFTReverseTransform(psFFT *fft        //!< FFT to apply
+                      );
+\end{verbatim}
+
+The data in Fourier space may be filtered using functions that return
+a multiplicative factor for a given position in Fourier space:
+
+\begin{verbatim}
+/** Apply filter function in fourier space */
+psFFT *
+psFFTFilter(psFFT *fft,                 //!< FFT to use (input and output)
+            float (*filterFunc)(int kx, int ky) //!< External filter function
+            );
+
+/** Apply complex filter function */
+psFFT *
+psFFTFilterComplex(psFFT *fft,          //!< FFT to use (input and output)
+                   float (*realFilterFunc)(int kx, int ky), //!< External filter function, real part
+                   float (*imagFilterFunc)(int kx, int ky) //!< External filter function, imaginary part
+                   );
+\end{verbatim}
+
+The cross-correlation function may be calculated from two Fourier transforms,
+and returns a new \code{psFFT}:
+
+\begin{verbatim}
+/** Calculate cross-correlation function */
+psFFT *
+psFFTCrossCorrelate(psFFT *out          //!< Output FFT (or NULL)
+                    psFFT *fft1, psFFT *fft2 //!< FFTs to use in cross-correlation
+                    );
+\end{verbatim}
+
+The power spectrum is calculated from the Fourier transform:
+
+\begin{verbatim}
+/** Calculate power spectrum */
+psFFT *
+psFFTPowerSpec(psFFT *fft               //!< FFT to use (input and output)
+               );
+\end{verbatim}
+
+In convolution, two Fourier transforms are multiplied.
+
+\begin{verbatim}
+/** Multiply two Fourier transforms, as for convolution */
+psFFT *
+psFFTMultiplyFT(psFFT *out,             //!< Output FFT (or NULL)
+                const psFFT *fft1, const psFFT *fft2 //!< FFTs to multiply
+                );
+\end{verbatim}
+
+And finally, the user may extract both the real data and the complex
+Fourier space data (without needing to destroy the \code{psFFT}).
+
+\begin{verbatim}
+/* Convert the real data in the FFT struct to an image again */
+psImage *
+psFFTGetImage(psImage *out,             //!< Image to write to (or NULL)
+              const psFFT *fft          //!< FFT to get image from
+              );
+
+/** Convert the Fourier transform data in the FFT struct to an image of complex numbers */
+psImage *
+psFFTGetFT(psImage *out,                //!< Image to write to (or NULL)
+           const psFFT *fft             //!< FFT to get Fourier transform from
+           );
 \end{verbatim}