I am currently working on a beam propagation code. The code basically consists of a transversal profile (a 2000x2000 array) which is propagated a small distance, then has its borders filtered, then propagated a small distance and so repeatedly. Here are two examples of filters I am using: I can provide code to the filters if needed, but the point is that they preserve the center of the images and absorb the borders of them; it is also important that they are continuous. Right now, what I am doing to apply the filters to an image U is: from numpy.fft import fft2, ifft2 Uf=fft2(U) #propagate Uf a small distance U=ifft2(Uf) U=U*F #F is the filter (a 2D array with the same dimensions as U) #repeat The Fourier transforms are the slowest part of my code, and I think I can do something more efficient by applying the filter in Fourier space. What is in my mind is the convolution theorem, which states that the transform of a product equals the convolution of the transforms (F[f.g]=F[f]*F[g]). So I thought that I could precompute the transform of my filter and then do the convolution between my image and the filter after each step, but it turns out that this approach is very slow. Is there a better way to do this? I read that doing two 1D convolutions instead of one 2D one can be faster, but I am not sure how to do it in my case. Continue reading...
