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. 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.