2D БПФ от MATLAB до Python - PullRequest
Новая
       29

2D БПФ от MATLAB до Python

2 голосов
/ 10 октября 2019

Большая часть моего опыта программирования на MATLAB, и я недавно начал знакомиться с Python.

Я натолкнулся на замечательный код MATLAB здесь , который относится к некоторым вещам, которые я хотел быработать, поэтому я попытался воссоздать его в 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. Любая идея, почемучастота y показана на 3 Гц, наоборот - 3 Гц
  2. Я не мог понять, что делал оригинальный комментатор в MATLAB с этими двумя строками: 
Nfft = 4 * 2 .^ nextpow2(size(im));
imF = fftshift(fft2(im, Nfft(1), Nfft(2))) / numel(im);

, чтовероятно, почему мой сигнал FFT такой интенсивный по отношению к фону. Мысли о том, как я могу настроить свое БПФ в Python?

1 Ответ

0 голосов
/ 10 октября 2019

У меня есть только частичный ответ.

Если вы внимательно посмотрите, цвета на синусоидальном изображении, сгенерированном с помощью кода Python, и цвета, сгенерированного с помощью кода Matlab, который вы связали, имеют инвертированный цвет (проверьте цвета полос ближе к краям и цвета нацветная полоса).

Это объясняет, почему у вас есть инвертированные цвета на графике FFT, и может быть, поэтому вы получили 3 Гц вместо -3 Гц. К сожалению, сейчас у меня нет доступа к компьютеру с Python, и я не смогу это проверить. Я думаю, это может быть хорошей идеей для начала поиска неисправностей.

РЕДАКТИРОВАТЬ:

Да, вы правы. Я полностью пропустил flipud в сценарии Matlab. Я не думаю, что ваш c расчет неверен. Самый простой способ проверить это - сохранить данные Matlab и импортировать их в Python.

В Matlab: 

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

Затем импортировать их в Python, используя scipy: 

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

Если последняя строка дает вам True, то это означает, что вычисления верны.

О графике, imshow предполагает, что источник (0-й индекс массива numpy)это верхний левый угол, который является нормой для изображений. Вы можете изменить это, используя ключевое слово origin='lower' в plt.imshow.

О fftshift, я думаю этот ответ на другой вопрос StackOverflow - это то, что вам нужно.

...