[python-dtcwt] 14/497: implement FFT-based 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 59a290e38243c6f2adbbbff0ba5e37bc4a39a46f
Author: Rich Wareham <rjw57 at cam.ac.uk>
Date:   Wed Aug 7 11:03:03 2013 +0100

    implement FFT-based convolution
---
 dtcwt/lowlevel.py | 40 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 39 insertions(+), 1 deletion(-)

diff --git a/dtcwt/lowlevel.py b/dtcwt/lowlevel.py
index c59c3f1..d52fa50 100644
--- a/dtcwt/lowlevel.py
+++ b/dtcwt/lowlevel.py
@@ -9,8 +9,46 @@ def _to_vertical_vector(v):
     else:
         return v
 
+def _centered(arr, newsize):
+    # Return the center newsize portion of the array.
+    # (Shamelessly cribbed from scipy.)
+    newsize = np.asarray(newsize)
+    currsize = np.array(arr.shape)
+    startind = (currsize - newsize) // 2
+    endind = startind + newsize
+    myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
+    return arr[tuple(myslice)]
+
 def _column_convolve(X, h):
-    return convolve2d(X, _to_vertical_vector(h), 'valid')
+    """Convolve the columns of X with h returning only the 'valid' section,
+    i.e. those values unaffected by zero padding.
+
+    """
+    h = h.flatten()
+    h_size = h.shape[0]
+    full_size = X.shape[0] + h_size - 1
+
+    # For small arrays, convolving directly is often faster
+    if full_size < 32:
+        return convolve2d(X, _to_vertical_vector(h), 'valid')
+
+    # Always use 2**n-sized FFT
+    fsize = 2 ** np.ceil(np.log2(full_size)).astype(int)
+
+    # Take FFT down columns
+    Xfft = np.fft.rfft(X, n=fsize, axis=0)
+
+    # Take FFT of input vector
+    hfft = np.fft.rfft(h, n=fsize)
+
+    # Column-wise multiply. I.e. scale rows of Xfft by hfft
+    Xfft = Xfft * hfft[:,np.newaxis]
+
+    # Invert
+    Xconv = np.fft.irfft(Xfft, n=fsize, axis=0)[:full_size,:].real
+    Xvalid = _centered(Xconv, (abs(X.shape[0] - h_size) + 1, X.shape[1]))
+
+    return Xvalid
 
 def reflect(x, minx, maxx):
     """Reflect the values in matrix x about the scalar values minx and maxx.

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