Is there a way to chose the x/y output axes range from np.fft2 ? I have a piece of code computing the diffraction pattern of an aperture. The aperture is defined in a 2k x 2k pixel array. The diffraction pattern is basically the inner part of the 2D FT of the aperture. The np.fft2 gives me an output array same size of the input but with some preset range of the x/y axes. Of course I can zoom in by using the image viewer, but I have already lost detail. What is the solution?
Thanks, Gert
import numpy as np
import matplotlib.pyplot as plt
r= 500
s= 1000
y,x = np.ogrid[-s:s+1, -s:s+1]
mask = x*x + y*y <= r*r
aperture = np.ones((2*s+1, 2*s+1))
aperture[mask] = 0
plt.imshow(aperture)
plt.show()
ffta= np.fft.fft2(aperture)
plt.imshow(np.log(np.abs(np.fft.fftshift(ffta))**2))
plt.show()
Unfortunately, much of the speed and accuracy of the FFT come from the outputs being the same size as the input.
The conventional way to increase the apparent resolution in the output Fourier domain is by zero-padding the input:
np.fft.fft2(aperture, [4 * (2*s+1), 4 * (2*s+1)])
tells the FFT to pad your input to be4 * (2*s+1)
pixels tall and wide, i.e., make the input four times larger (sixteen times the number of pixels).Begin aside I say "apparent" resolution because the actual amount of data you have hasn't increased, but the Fourier transform will appear smoother because zero-padding in the input domain causes the Fourier transform to interpolate the output. In the example above, any feature that could be seen with one pixel will be shown with four pixels. Just to make this fully concrete, this example shows that every fourth pixel of the zero-padded FFT is numerically the same as every pixel of the original unpadded FFT:
(That
relerr
is just a simple relative error, which you want to be close to machine precision, around2e-16
. The largest error between every 4th sample of the zero-padded FFT and the unpadded FFT is6e-12
which is quite close to machine precision, meaning these two arrays are nearly numerically equivalent.) End asideZero-padding is the most straightforward way around your problem. But it does cost you a lot of memory. And it is frustrating because you might only care about a tiny, tiny part of the transform. There's an algorithm called the chirp z-transform (CZT, or colloquially the "zoom FFT") which can do this. If your input is
N
(for you2*s+1
) and you want justM
samples of the FFT's output evaluated anywhere, it will compute three Fourier transforms of sizeN + M - 1
to obtain the desiredM
samples of the output. This would solve your problem too, since you can ask forM
samples in the region of interest, and it wouldn't require prohibitively-much memory, though it would need at least 3x more CPU time. The downside is that a solid implementation of CZT isn't in Numpy/Scipy yet: see the scipy issue and the code it references. Matlab's CZT seems reliable, if that's an option; Octave-forge has one too and the Octave people usually try hard to match/exceed Matlab.But if you have the memory, zero-padding the input is the way to go.