[python-dtcwt] 79/497: move from FFT based method for convolution to a direct convolution

[Ghislain Vaillant]

This is an automated email from the git hooks/post-receive script.

ghisvail-guest pushed a commit to branch debian/sid
in repository python-dtcwt.

commit 6436b41dd5952d639f8baf8aa8ad931cca92fc41
Author: Rich Wareham <rjw57 at cam.ac.uk>
Date:   Thu Aug 8 21:56:18 2013 +0100

    move from FFT based method for convolution to a direct convolution
    
    Since h is likely to be so small, it makes sense to do a direct
    convolution.
---
 dtcwt/lowlevel.py | 34 ++++++++++++----------------------
 1 file changed, 12 insertions(+), 22 deletions(-)

diff --git a/dtcwt/lowlevel.py b/dtcwt/lowlevel.py
index 008bf5f..4e94c43 100644
--- a/dtcwt/lowlevel.py
+++ b/dtcwt/lowlevel.py
@@ -1,4 +1,5 @@
 import numpy as np
+from six.moves import xrange
 
 def as_column_vector(v):
     """Return *v* as a column vector with shape (N,1).
@@ -28,34 +29,23 @@ def _column_convolve(X, h):
     """Convolve the columns of *X* with *h* returning only the 'valid' section,
     i.e. those values unaffected by zero padding.
 
-    """
-
-    # This function should give the same result as:
-    #
-    # from scipy.signal import convolve2d
-    # return convolve2d(X, as_column_vector(h), 'valid')
+    We assume that h is small and so direct convolution is the most efficient.
 
+    """
+    Xshape = np.asarray(X.shape)
     h = h.flatten()
     h_size = h.shape[0]
-    full_size = X.shape[0] + h_size - 1
 
-    # Always use 2**n-sized FFT
-    fsize = 2 ** np.ceil(np.log2(full_size)).astype(int)
-
-    # Take FFT down columns
-    Xfft = _rfft(X, n=fsize, axis=0)
-
-    # Take FFT of input vector
-    hfft = _rfft(h, n=fsize, axis=0)
-
-    # Column-wise multiply. I.e. scale rows of Xfft by hfft
-    Xfft = Xfft * hfft[:,np.newaxis]
+    full_size = X.shape[0] + h_size - 1
+    Xshape[0] = full_size
 
-    # Invert
-    Xconv = _irfft(Xfft, n=fsize, axis=0)[:full_size,:].real
-    Xvalid = _centered(Xconv, (abs(X.shape[0] - h_size) + 1, X.shape[1]))
+    out = np.zeros(Xshape, dtype=X.dtype)
+    for idx in xrange(h_size):
+        out[idx:(idx+X.shape[0]),...] += X * h[idx]
 
-    return Xvalid
+    outShape = Xshape.copy()
+    outShape[0] = abs(X.shape[0] - h_size) + 1
+    return _centered(out, outShape)
 
 def reflect(x, minx, maxx):
     """Reflect the values in matrix *x* about the scalar values *minx* and

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/python-dtcwt.git