2D FFT from MATLAB to Python

user_12192138 0
Viewed 641

Most of my programming experience is in MATLAB and I recently started get familiar with Python.

I came across some great MATLAB code here that pertains to some things I'd like to work with, so I've tried to recreate it in Python:

import numpy as np
import math
import matplotlib.pyplot as plt

x = np.linspace(-2, 2, 100) # seconds
y = np.linspace(-3, 3, 200) # seconds

xFreq = 2; # Hz
yFreq = -3; # Hz

a = np.matrix(np.matrix(np.exp(2j * np.pi * y * yFreq)))
b = np.matrix(np.exp(2j * np.pi * np.matrix(x).T * xFreq))

c = np.dot(b,a).T

plt.imshow(c.real, cmap='gray', extent = [min(x), max(x), min(y), max(y)], aspect=2/3);
plt.colorbar()
plt.xlabel('x (Sec)')
plt.ylabel('y (Sec)')
plt.show()

2D sinuosidal

nfftx = len(x);
fs = 1/np.diff(x)[0]; 
fx = np.linspace(-1,1,nfftx) * fs/2;

nffty = len(y);
fs = 1/np.diff(y)[0]; 
fy = np.linspace(-1,1,nffty) * fs/2;

imF = np.fft.fftshift(np.fft.fft2(c))/np.size(c)

plt.figure()
plt.title("FFT (real)")
plt.imshow(np.real(imF), cmap='gray', extent = [min(fx), max(fx), min(fy), max(fy)], aspect=2/3)
plt.xlabel('fx (Hz)')
plt.ylabel('fy (Hz)')

FFT of the image

  1. Any idea why the y frequency is shown at 3 Hz vice -3 Hz
  2. I couldn't understand what the original commentator was doing in MATLAB with these two lines:
Nfft = 4 * 2 .^ nextpow2(size(im));
imF = fftshift(fft2(im, Nfft(1), Nfft(2))) / numel(im);

which is likely why my FFT signal is so intense relative to background. Thoughts on how I could adjust my FFT in Python?

fft image-processing matlab python
1 Answers

I only have a partial answer.

If you look closely, the colors on the sinusoidal image generated with your Python code and the one generated with Matlab code you linked have inverted color (check the colors of the stripes closer to edges, and the colors on the color bar).

That explains why you have inverted colors on the FFT plot, and may be why you got 3 Hz, instead of -3 Hz. Unfortunately, I do not have access to a computer with Python right now and won't be able to verify this. I guess this may be a good thing to start troubleshooting with.

EDIT:

Yes, you are right. I completely missed the flipud in the Matlab script. I do no think your c calculation is wrong. The easiest way to check that is to save the Matlab data and import it to Python.

In Matlab:

save('data.mat', 'im')

Then import that to Python using scipy:

im_matlab = scipy.io.loadmat('data.mat')['im']
np.all(np.isclose(im_matlab, im))

If the last line gives you True, then that means the calculations are correct.

About the plot, imshow assumes that the origin (0-th index of the numpy array) is the top left corner, which is the norm for images. You can change this by using origin='lower' keyword in plt.imshow.

About fftshift, I think this answer to a different StackOverflow question is what you need.

user_83716910
Related