Added fft() and ifft()

dlib/matrix.h

@@ -12,6 +12,7 @@
 #include "matrix/symmetric_matrix_cache.h"
 #include "matrix/matrix_conv.h"
 #include "matrix/matrix_read_from_istream.h"
+#include "matrix/matrix_fft.h"
 
 #ifdef DLIB_USE_BLAS
 #include "matrix/matrix_blas_bindings.h"

dlib/matrix/matrix_fft.h

+// Copyright (C) 2013  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_FFt_H__
+#define DLIB_FFt_H__
+
+#include "matrix_fft_abstract.h"
+#include "matrix_utilities.h"
+#include "../hash.h"
+#include "../algs.h"
+
+#ifdef DLIB_USE_FFTW
+#include <fftw3.h>
+#endif // DLIB_USE_FFTW
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    namespace impl
+    {
+        inline unsigned long reverse_bits (
+            unsigned long val,
+            unsigned long num
+        )
+        {
+            unsigned long temp = 0;
+            for (unsigned long i = 0; i < num; ++i)
+            {
+                temp <<= 1;
+                temp |= val&0x1;
+                val >>= 1;
+            }
+            return temp;
+        }
+
+        template <typename T, long NR, long NC, typename MM, typename L>
+        void permute (
+            const matrix<std::complex<T>,NR,NC,MM,L>& data, 
+            matrix<std::complex<T>,NR,NC,MM,L>& outdata 
+        )  
+        {
+            outdata.set_size(data.size());
+            if (data.size() == 0)
+                return;
+
+            const unsigned long num = static_cast<unsigned long>(std::log(data.size())/std::log(2.0) + 0.5);
+            for (unsigned long i = 0; i < (unsigned long)data.size(); ++i)
+            {
+                outdata(impl::reverse_bits(i,num)) = data(i);
+            }
+        }
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    inline bool is_power_of_two (
+        const unsigned long& value
+    )
+    {
+        if (value == 0)
+            return true;
+        else
+            return count_bits(value) == 1;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename T, long NR, long NC, typename MM, typename L>
+    matrix<std::complex<T>,NR,NC,MM,L> fft (
+        const matrix<std::complex<T>,NR,NC,MM,L>& data
+    )  
+    {
+        if (data.size() == 0)
+            return data;
+
+        // make sure requires clause is not broken
+        DLIB_CASSERT(is_vector(data) && is_power_of_two(data.size()),
+            "\t void ifft(data)"
+            << "\n\t data must be a vector with a size that is a power of two."
+            << "\n\t is_vector(data): " << is_vector(data)
+            << "\n\t data.size():     " << data.size()
+            );
+
+        matrix<std::complex<T>,NR,NC,MM,L> outdata(data);
+
+        const long half = outdata.size()/2;
+
+        typedef std::complex<T> ct;
+        matrix<ct,0,1,MM,L> twiddle_factors(half);
+
+        // compute the complex root of unity w
+        const T temp = -2.0*pi/outdata.size();
+        ct w = ct(std::cos(temp),std::sin(temp));
+
+        ct w_pow = 1;
+
+        // compute the twiddle factors
+        for (long j = 0; j < twiddle_factors.size(); ++j)
+        {
+            twiddle_factors(j) = w_pow; 
+            w_pow *= w;
+        }
+
+
+        // now compute the decimation in frequency.  This first
+        // outer loop loops log2(outdata.size()) number of times
+        long skip = 1;
+        for (long step = half; step != 0; step >>= 1)
+        {
+            // do blocks of butterflies in this loop
+            for (long j = 0; j < outdata.size(); j += step*2)
+            {
+                // do step butterflies
+                for (long k = 0; k < step; ++k)
+                {
+                    const long a_idx = j+k;
+                    const long b_idx = j+k+step;
+                    const ct a = outdata(a_idx) + outdata(b_idx);
+                    const ct b = (outdata(a_idx) - outdata(b_idx))*twiddle_factors(k*skip);
+                    outdata(a_idx) = a;
+                    outdata(b_idx) = b;
+                }
+            }
+            skip *= 2;
+        }
+
+        matrix<std::complex<T>,NR,NC,MM,L> outperm;
+        impl::permute(outdata, outperm);
+        return outperm;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename T, long NR, long NC, typename MM, typename L>
+    matrix<std::complex<T>,NR,NC,MM,L> ifft (
+        const matrix<std::complex<T>,NR,NC,MM,L>& data 
+    )  
+    {
+        if (data.size() == 0)
+            return data;
+
+        // make sure requires clause is not broken
+        DLIB_CASSERT(is_vector(data) && is_power_of_two(data.size()),
+            "\t void ifft(data)"
+            << "\n\t data must be a vector with a size that is a power of two."
+            << "\n\t is_vector(data): " << is_vector(data)
+            << "\n\t data.size():     " << data.size()
+            );
+
+
+        matrix<std::complex<T>,NR,NC,MM,L> outdata;
+        impl::permute(data,outdata);
+
+        const long half = outdata.size()/2;
+
+        typedef std::complex<T> ct;
+        matrix<ct,0,1,MM,L> twiddle_factors(half);
+
+        // compute the complex root of unity w
+        const T temp = 2.0*pi/outdata.size();
+        ct w = ct(std::cos(temp),std::sin(temp));
+
+        ct w_pow = 1;
+
+        // compute the twiddle factors
+        for (long j = 0; j < twiddle_factors.size(); ++j)
+        {
+            twiddle_factors(j) = w_pow; 
+            w_pow *= w;
+        }
+
+        // now compute the inverse decimation in frequency.  This first
+        // outer loop loops log2(outdata.size()) number of times
+        long skip = half;
+        for (long step = 1; step <= half; step <<= 1)
+        {
+            // do blocks of butterflies in this loop
+            for (long j = 0; j < outdata.size(); j += step*2)
+            {
+                // do step butterflies
+                for (long k = 0; k < step; ++k)
+                {
+                    const long a_idx = j+k;
+                    const long b_idx = j+k+step;
+                    outdata(b_idx) *= twiddle_factors(k*skip);
+                    const ct a = outdata(a_idx) + outdata(b_idx);
+                    const ct b = outdata(a_idx) - outdata(b_idx);
+                    outdata(a_idx) = a;
+                    outdata(b_idx) = b;
+                }
+            }
+            skip /= 2;
+        }
+
+        outdata /= outdata.size();
+        return outdata;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+#ifdef DLIB_USE_FFTW
+
+    inline matrix<std::complex<double>,0,1> fft(
+        const matrix<std::complex<double>,0,1>& data
+    )
+    {
+        // make sure requires clause is not broken
+        DLIB_CASSERT(is_vector(data) && is_power_of_two(data.size()),
+            "\t void fft(data)"
+            << "\n\t data must be a vector with a size that is a power of two."
+            << "\n\t is_vector(data): " << is_vector(data)
+            << "\n\t data.size():     " << data.size()
+            );
+
+        matrix<std::complex<double>,0,1> m2(data.size());
+        fftw_complex *in, *out;
+        fftw_plan p;
+        in = (fftw_complex*)&data(0);
+        out = (fftw_complex*)&m2(0);
+        p = fftw_plan_dft_1d(data.size(), in, out, FFTW_FORWARD, FFTW_ESTIMATE);
+        fftw_execute(p); 
+        fftw_destroy_plan(p);
+        return m2;
+    }
+
+    inline matrix<std::complex<double>,0,1> ifft(
+        const matrix<std::complex<double>,0,1>& data
+    )
+    {
+        // make sure requires clause is not broken
+        DLIB_CASSERT(is_vector(data) && is_power_of_two(data.size()),
+            "\t void ifft(data)"
+            << "\n\t data must be a vector with a size that is a power of two."
+            << "\n\t is_vector(data): " << is_vector(data)
+            << "\n\t data.size():     " << data.size()
+            );
+
+        matrix<std::complex<double>,0,1> m2(data.size());
+        fftw_complex *in, *out;
+        fftw_plan p;
+        in = (fftw_complex*)&data(0);
+        out = (fftw_complex*)&m2(0);
+        p = fftw_plan_dft_1d(data.size(), in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
+        fftw_execute(p); 
+        fftw_destroy_plan(p);
+        return m2/data.size();
+    }
+
+#endif // DLIB_USE_FFTW
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_FFt_H__
+

dlib/matrix/matrix_fft_abstract.h

+// Copyright (C) 2013  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_FFt_ABSTRACT_H__
+#ifdef DLIB_FFt_ABSTRACT_H__
+
+#include "matrix_abstract.h"
+#include "../algs.h"
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    bool is_power_of_two (
+        const unsigned long& value
+    );
+    /*!
+        ensures
+            - returns true if value contains a power of two and false otherwise.  As a
+              special case, we also consider 0 to be a power of two.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T, 
+        long NR, 
+        long NC, 
+        typename MM, 
+        typename L
+        >
+    matrix<std::complex<T>,NR,NC,MM,L> fft (
+        const matrix<std::complex<T>,NR,NC,MM,L>& data
+    );
+    /*!
+        requires
+            - is_vector(data) == true
+            - is_power_of_two(data.size()) == true
+        ensures
+            - Computes the discrete Fourier transform of the given data vector and
+              returns it.  In particular, we return a matrix D such that:
+                - D.nr() == data.nr()
+                - D.nc() == data.nc()
+                - D(0) == the DC term of the Fourier transform.
+                - starting with D(0), D contains progressively higher frequency components
+                  of the input data.
+                - ifft(D) == D
+            - if DLIB_USE_FFTW is #defined then this function will use the very fast fftw
+              library when given double precision matrices instead of dlib's default fft
+              implementation.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T, 
+        long NR, 
+        long NC, 
+        typename MM, 
+        typename L
+        >
+    matrix<std::complex<T>,NR,NC,MM,L> ifft (
+        const matrix<std::complex<T>,NR,NC,MM,L>& data 
+    );
+    /*!
+        requires
+            - is_vector(data) == true
+            - is_power_of_two(data.size()) == true
+        ensures
+            - Computes the inverse discrete Fourier transform of the given data vector and
+              returns it.  In particular, we return a matrix D such that:
+                - D.nr() == data.nr()
+                - D.nc() == data.nc()
+                - fft(D) == data 
+            - if DLIB_USE_FFTW is #defined then this function will use the very fast fftw
+              library when given double precision matrices instead of dlib's default fft
+              implementation.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_FFt_ABSTRACT_H__
+

dlib/test/CMakeLists.txt

@@ -44,6 +44,7 @@ set (tests
    entropy_coder.cpp
    entropy_encoder_model.cpp
    example_args.cpp
+   fft.cpp
    filtering.cpp
    find_max_factor_graph_nmplp.cpp
    find_max_factor_graph_viterbi.cpp

dlib/test/makefile

@@ -59,6 +59,7 @@ SRC += ekm_and_lisf.cpp
 SRC += empirical_kernel_map.cpp
 SRC += entropy_coder.cpp
 SRC += entropy_encoder_model.cpp
+SRC += fft.cpp
 SRC += filtering.cpp
 SRC += find_max_factor_graph_nmplp.cpp
 SRC += find_max_factor_graph_viterbi.cpp