Расчет лапласиана в реальном pyFFTW - PullRequest
Новая
       42

Расчет лапласиана в реальном pyFFTW

3 голосов
/ 28 июня 2019

Для прямого (многомерного) алгоритма FFTW вы можете указать, что вход numpy.ndarray является действительным, а вывод должен быть сложным.Это делается при создании выровненных по байту массивов, которые идут в аргументах fft_object: 

import numpy as np
import pyfftw

N = 256  # Input array size (preferrably 2^{a}*3^{b}*5^{c}*7^{d}*11^{e}*13^{f}, (e+f = 0,1))
dx = 0.1  # Spacing between mesh points
a = pyfftw.empty_aligned((N, N), dtype='float64')
b = pyfftw.empty_aligned((N, N//2+1), dtype='complex128')
fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_FORWARD')

Выходной массив не симметричен, а вторая ось усекается до положительных частот.Для сложного БПФ вы можете вычислить лапласиан со следующим np.ndarray

kx, ky = np.meshgrid(np.fft.fftfreq(N, dx), np.fft.fftfreq(N, dx))  # Wave vector components
k2 = -4*np.pi**2*(kx*kx+ky*ky)  # np.ndarray for the Laplacian operator in "frequency space"

Как это будет сделано в усеченном случае?Я думал об использовании: 

kx, ky = np.meshgrid(np.fft.fftfreq(N//2+1, dx), np.fft.fftfreq(N, dx))  # The axes conven-
#                                                                        tions are different

Но будет ли это действительно работать?Похоже, что он пренебрегает отрицательными частотами в направлении "у".

1 Ответ

1 голос
/ 28 июня 2019

Я не знаком с pyfftw, но с модулем numpy.fft он будет работать просто отлично (при условии, что вы используете rfftfreq, как указано в комментариях).

Напомним: длявещественный массив, a, преобразование Фурье, b, обладает свойством Эрмита: b(-kx,-ky) - комплексное сопряжение b(kx,ky).Реальная версия forward fft отбрасывает (большую часть) избыточную информацию, пропуская отрицательные ky s.Реальная версия обратного БПФ предполагает, что значения на пропущенных частотах можно найти путем комплексного сопряжения соответствующих элементов.

Если бы вы использовали комплексное БПФ и сохранили все частоты, -k2 * b все равно будет иметьЭрмитово-подобная собственность.Таким образом, предположение, сделанное реальным обратным fft, все еще остается верным и даст правильный ответ.

Я думаю, с pyfftw все будет работать нормально, если вы укажете массив float64 правильного размера длявыход для случая direction=FFT_BACKWARD.

...