2

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?

NewtoPy
  • 23
  • 3
  • The plotting with imshow is upside down (use origin='lower' for your case). Your data are synthetic, so you really don't expect anything but a white dot. Add some noise to your model. Also you can plot a root of the spectrum or the log of it, when you expect but can't see finer structure. – roadrunner66 Oct 12 '19 at 00:34
  • This was it. Thank you! – NewtoPy Oct 15 '19 at 11:30

1 Answers1

-1

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.

Kani
  • 469
  • 4
  • 13
  • I’m almost positive that the discrepancy in color you see is due to the OP inverting the color scale on his/her MATLAB script. I should say that I suspect the discrepancy in my y Hz is likely due to the matrix multiplication happening in creating the variable “c”. – NewtoPy Oct 10 '19 at 10:10
  • Yes, you are right. I completely missed `flipud` in the OP. I edited my answer. – Kani Oct 10 '19 at 13:16
  • Unfortunately, I don’t have a MATLAB license anymore. – NewtoPy Oct 10 '19 at 13:37
  • It seems that origin = 'lower' was the fix! I verified that I had calculated the correct frequencies by finding the indices of the peak signal and finding the corresponding fx and fy. – NewtoPy Oct 15 '19 at 11:27
  • I finally got access to a computer that has both Matlab and Python installed today. The Python calculation of `c` exactly matches the Matlab one too. – Kani Oct 15 '19 at 14:09