1
$\begingroup$

I've been working on a small project, and as a part of this I need a function that can compute the structure function of the input data. In the process of writing this part/testing it I found the results were actually a poor match to what I would expect based on an analytic solution. Upon some further investigation, I've found this mismatch applies even to 2D convolutions and I have been completely stuck trying to debug this.

The test case that I am running is two rectangular functions convolved with one another, my expected result is a 2D triangle function (a pyramid essentially). But the difference between the analytic output and the FFT convolution output is large enough that I am concerned something is wrong. Furthermore, when I compute a direct convolution of the two inputs using scipy.convolve I get yet another output, different from the other two in it's own unique way!

To demonstrate the differences I'm seeing, I've plotted: my FFT Convolution's output, the analytic solution, the difference between them, and the difference between scipy.convolve's output and the analytic solution.Comparison of the three different convolution outputs My question then boils down to a couple things:

  • How large of a difference is expected between the FFT convolution and the analytic solution?
    • I would have expected something small (of order 1e-15 or less) but the peak difference I am seeing is 0.12. It also follows a very particular pattern, the vast majority of the output is close enough (difference of ~1.0e-16).
  • Is there a key step or concept that I am missing here?
    • I've been taking the FFT of each input, multiplying them together, and then taking the inverse FFT (with some scaling based on pixel size).

The code that produces the figure can be found in this relatively simple jupyter notebook on my github.

$\endgroup$
2
$\begingroup$

After playing with the notebook for a while, I figured out the following:

  1. The result from scipy is shifted by one pixel in each axis (probably due to its definition of the "same" mode). Or, alternatively, both the analytic result and pyfftw results are shifted the opposite way.

    After working around that, the result is exactly (within error) the same as the pyfftw result. I did this as follows (though this drops some data on the edges):

B_padded = np.pad(B, ((0, 1), (0, 1)), 'edge')[1:, 1:]
C1 = scipy.signal.convolve(A, B_padded, 'same', 'direct') * delta**2
  1. The analytic solution (which peaks at 4.0) is not quite right for your input data (31x31 pixels at 1.0, edges at 0.5, and corners at 0.25). The peak value (where the image and the window align perfectly) should be:
box_bulk = 1.0**2 * (31 * 31)
box_edges = 0.5**2 * 31 * 4
box_corners = 0.25**2 * 4
peak = (box_bulk + box_edges + box_corners) / 1024 * 4  # ~= 3.876

which is the answer produced by both scipy and pyfftw.

Gist (also includes plotting tweaks).

(disclaimer: friend of asker)

| improve this answer | |
$\endgroup$
  • $\begingroup$ That's great, the scipy convolve agreeing with the FFT convolve seems like a great step forward. I still am curious what the analytic solution should be though, since changing its overall normalization just shifts the disagreement from being along the cross to in the quadrants. $\endgroup$ – LocalMin Mar 9 at 3:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.