Python: как интерполировать «неструктурированные» данные двумерного преобразования Фурье - PullRequest
Новая
       33

Python: как интерполировать «неструктурированные» данные двумерного преобразования Фурье

0 голосов
/ 06 июля 2019

Моя цель - интерполировать дискретное непрерывное двумерное преобразование Фурье функции.Кажется, проблема в том, что частоты в каждом измерении не выводятся в строго возрастающем порядке (см. здесь ).

Функция fft.fft2 принимает двумерный массив, где в моем случаемассив (назовем его A) структурирован так, что A[i][j] = fun(x[i], y[j]), fun - функция, подлежащая преобразованию.После применения fft.fft2 к A на выходе получается массив F того же размера, что и исходный массив, так что частотная координата, соответствующая F[i][j], равна (w_x[i], w_y[j]), где w_x = fft.fftfreq(F.shape[0]) и w_y = fft.fftfreq(F.shape[1]),оба из них являются одномерными массивами, которые не в порядке возрастания.

За wx и wy Я хочу интерполировать F (скажем, функции finterp) так, чтобы интерполированное значение быловозвращается при вызове finterp(w_x, w_y), w_x и w_y в пределах домена wx и диапазона wy, но в противном случае произвольно.Я рассмотрел варианты интерполяции, доступные через scipy.interpolate , но мне не кажется, что кто-либо из них может иметь дело с этим типом структуры данных (координатные оси определяются как внешние1D массивы порядка и значения функций, находящиеся в 2D массиве).

Это немного абстрактно, поэтому здесь я составил простой пример, похожий по структуре на приведенный выше.Предположим, что мы хотим построить непрерывную функцию f(x, y) = x + y в области x = [-1, 1] и y = [-1, 1], учитывая следующие данные: 

import numpy as np

# note that below z[i][j] corresponds to what we want f(x[i], y[j]) to be

x = np.array([0, 1, -1])
y = np.array([0, 1, -1])
z = np.array([0, 1, -1],[1, 2, 0],[-1, 0, -2])

z[i][j], как мы знаем, соответствует функции, оцененной в x[i], y[j],Как можно (а) напрямую интерполировать эти данные, учитывая их исходную структуру, или (б) переставить данные так, чтобы x и y были в порядке возрастания, а упорядоченное z было таким, чтобы z[i][j]равно функции, оцененной при перестановке x[i], y[j]?

1 Ответ

1 голос
/ 08 июля 2019

Следующий код показывает, как использовать fftshift для изменения выходных данных fft2 и fftfreq, чтобы оси частот монотонно увеличивались.После применения fftshift вы можете использовать массивы для интерполяции.Я добавил отображение массивов, чтобы вы могли убедиться, что сами данные не изменились.Начало координат смещено из верхнего левого угла в середину массива, перемещая отрицательные частоты с правой стороны на левую. 

import numpy as np
import matplotlib.pyplot as pp

x = np.array([0, 1, -1])
y = np.array([0, 1, -1])
z = np.array([[0, 1, -1],[1, 2, 0],[-1, 0, -2]])
f = np.fft.fft2(z)
w_x = np.fft.fftfreq(f.shape[0])
w_y = np.fft.fftfreq(f.shape[1])

pp.figure()
pp.imshow(np.abs(f))
pp.xticks(np.arange(0,len(w_x)), np.round(w_x,2))
pp.yticks(np.arange(0,len(w_y)), np.round(w_y,2))

f = np.fft.fftshift(f)
w_x = np.fft.fftshift(w_x)
w_y = np.fft.fftshift(w_y)

pp.figure()
pp.imshow(np.abs(f))
pp.xticks(np.arange(0,len(w_x)), np.round(w_x,2))
pp.yticks(np.arange(0,len(w_y)), np.round(w_y,2))
pp.show()

Альтернативный подход состоит в том, чтобы не использовать fftfreq, чтобы определить свои частоты, но вычислить их вручную.БПФ по умолчанию вычисляет ДПФ для k=[0..N-1].Из-за периодичности, когда DFT на k равен DFT на k+N и k-N, его выход часто интерпретируется как k=[N//2...(N-1)//2] (но по-разному соответствует k=[0..N-1]);это k, которое fftfreq возвращает (возвращает k/N).

Таким образом, вместо этого вы можете сказать 

N = f.shape[0]
w_x = np.linspace(0, N, N, endpoint=False) / N

Теперь у вас нет отрицательных частот,и вместо этого иметь частоты в диапазоне [0,N-1]/N.

...