Added xcorr_fft()

dlib/matrix/matrix_conv.h

@@ -5,6 +5,7 @@
 
 #include "matrix_conv_abstract.h"
 #include "matrix.h"
+#include "matrix_fft.h"
 
 namespace dlib
 {
@@ -112,6 +113,55 @@ namespace dlib
         return matrix_op<op>(op(m1.ref(),m2.ref()));
     }
 
+// ----------------------------------------------------------------------------------------
+
+    namespace impl
+    {
+        inline size_t bounding_power_of_two (
+            size_t n
+        )
+        {
+            size_t s = 1;
+            for (unsigned int i = 0; i < sizeof(s)*8 && s < n; ++i)
+                s <<= 1;
+            return s;
+        }
+    }
+
+    template <
+        typename EXP1, 
+        typename EXP2
+        >
+    typename EXP1::matrix_type xcorr_fft(
+        const matrix_exp<EXP1>& u,
+        const matrix_exp<EXP2>& v
+    )
+    {
+        COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type, typename EXP2::type>::value == true));
+        using T = typename EXP1::type;
+        COMPILE_TIME_ASSERT((is_same_type<double,T>::value || is_same_type<float,T>::value || is_same_type<long double,T>::value ));
+
+        const long pad_nr = impl::bounding_power_of_two(u.nr() + v.nr() - 1);
+        const long pad_nc = impl::bounding_power_of_two(u.nc() + v.nc() - 1);
+
+        matrix<std::complex<T>> U(pad_nr, pad_nc), V(pad_nr,pad_nc);
+
+        U = 0;
+        V = 0;
+        set_subm(U,U.nr()-u.nr(),U.nc()-u.nc(),u.nr(),u.nc()) = u;
+        set_subm(V,get_rect(v)) = v;
+
+        fft_inplace(U);
+        fft_inplace(V);
+
+        return subm(real(ifft(pointwise_multiply(U, conj(V)))),
+            U.nr()-u.nr()-v.nr()+1, 
+            U.nc()-u.nc()-v.nc()+1, 
+            u.nr()+v.nr()-1,
+            u.nc()+v.nc()-1
+        );
+    }
+
 // ----------------------------------------------------------------------------------------
 // ----------------------------------------------------------------------------------------
 

dlib/matrix/matrix_conv_abstract.h

@@ -45,6 +45,23 @@ namespace dlib
                 - R.nc() == m1.nc()+m2.nc()-1
     !*/
 
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp xcorr_fft (
+        const matrix_exp& m1,
+        const matrix_exp& m2
+    );
+    /*!
+        requires
+            - m1 and m2 both contain elements of the same type
+            - m1 and m2 contain real or complex values and must be double, float, or long
+              double valued. (e.g. not integers)
+        ensures
+            - This function is identical to xcorr() except that it uses a fast Fourier
+              transform to do the convolution and is therefore much faster when both m1 and
+              m2 are large.
+    !*/
+
 // ----------------------------------------------------------------------------------------
 
     const matrix_exp conv_same (

dlib/test/matrix4.cpp

@@ -794,9 +794,9 @@ namespace
         DLIB_TEST(temp3 == temp);
 
 
+        dlib::rand rnd;
         for (int i = 0; i < 3; ++i)
         {
-            dlib::rand rnd;
             matrix<complex<int> > a, b;
             a = complex_matrix(matrix_cast<int>(round(20*randm(2,7,rnd))), 
                                matrix_cast<int>(round(20*randm(2,7,rnd))));
@@ -807,6 +807,20 @@ namespace
             DLIB_TEST(xcorr_valid(a,b) == conv_valid(a, flip(conj(b))));
             DLIB_TEST(xcorr_same(a,b)  == conv_same(a, flip(conj(b))));
         }
+
+
+        for (int i = 0; i < 30; ++i)
+        {
+            auto nr1 = rnd.get_integer_in_range(1,30);
+            auto nc1 = rnd.get_integer_in_range(1,30);
+            auto nr2 = rnd.get_integer_in_range(1,30);
+            auto nc2 = rnd.get_integer_in_range(1,30);
+            matrix<double> a, b;
+            a = randm(nr1,nc1,rnd);
+            b = randm(nr2,nc2,rnd);
+
+            DLIB_TEST(max(abs(xcorr(a,b) - xcorr_fft(a,b))) < 1e-12);
+        }
     }
 
     void test_complex()