CODEWITHSUNDEEP
pythonmatlabfft
user3090569
user3090569

Reputation: 21

Converting a Matlab fft2 diffraction example into Python

I trying to teach myself some FFT basics using python. At the moment I'm trying to reproduce a matlab FFT2 diffraction result using python and associated libraries.

The example I'm trying to simulate can be found here: http://www.mathworks.ch/ch/help/matlab/math/fast-fourier-transform-fft.html#brfb2vw-3

At the bottom of the page you will see this matlab code:

D1 = fft2(M);
D2 = fftshift(D1);

imagesc(abs(D2))
axis image
colormap(hot)
title('{\bf Diffraction Pattern}')

This turns converts a circular aperture image (they actually build this as an array if it makes a difference) into a diffraction pattern.

[Circular aperture] www.mathworks.ch/help/releases/R2013b/matlab/math/aperture.gif

[Matlab calculation result] www.mathworks.ch/help/releases/R2013b/matlab/math/diffraction1.gif

The problem is that when I take their circular aperture image (as a gif graphic) and apply the equivalent scipy FFT2 instructions I get a different result. Why this difference and what is the correct way to create the diffraction pattern. Are there similar gotchas for going back from the diffraction pattern back to the image that made it? Does this only work if one builds the aperture as an array like this?

My python code is as follows:

import sys
import numpy as np
import pylab as py
from scipy import misc, fftpack 

image = misc.imread(sys.argv[1])

D1 = fftpack.fft2(image)
D2 = fftpack.fftshift(D1)

abs_image = np.abs(D2)
py.imshow(abs_image)
py.show()

I downloaded and used their circular aperture image as input (after cropping of the scales etc), but I get a coloured cross hair on a blue background as a result. There is something that looks like a small dot in the middle of the cross hair. Is this result due to the fact that I am using the downloaded circular aperture image with its square edges or is it due to the code?

Many thanks in advance

Upvotes: 2

Views: 3286

Answers (1)

Stefan Seemayer
Stefan Seemayer

Reputation: 2067

My version of scipy.misc.imread will give you a W*H*4 array of RGBA colors, not intensity values and the FFT on it will probably give you different results.

I've ported the Matlab example code for generating M into numpy:

import numpy as np
import pylab as py
from scipy import misc, fftpack

n = 2**10
I = np.arange(1, n)
x = I - n / 2
y = n / 2 - I

R = 10

X = x[:, np.newaxis]
Y = y[np.newaxis, :]

M = X**2 + Y**2 < R**2

D1 = fftpack.fft2(M)
D2 = fftpack.fftshift(D1)

abs_image = np.abs(D2)
py.imshow(abs_image)
py.show()

This gives a better result:

Result

Upvotes: 2

Related Questions