Что не так с этим преобразованием Фурье?(в питоне) - PullRequest
Новая
       48

Что не так с этим преобразованием Фурье?(в питоне)

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

Как могли бы сказать люди, которые читали мои предыдущие посты на этом сайте, я пытался реализовать решатель PDE, который использует FFT в Python. Программная часть в основном проработана, но программа выдает (очень подходящую для этого сайта) ошибку переполнения (в основном она очень сильно увеличивается, пока не станет NaN).

После исключения всех других возможностей я связал проблему с БПФ и способом, которым я пытаюсь сделать производные, поэтому я решил протестировать два разных БПФ (numpy 'fft модуль и pyFFTW пакет) со следующим кодом: 

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


def fftw_(y: np.ndarray) -> np.ndarray:
    a = pyfftw.empty_aligned((N, N), dtype=np.float64)
    b = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
    fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_FORWARD', flags=('FFTW_MEASURE',), threads=12)
    y_hat = fft_object(y)
    return y_hat


def ifftw_(y_hat: np.ndarray) -> np.ndarray:
    a = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
    b = pyfftw.empty_aligned((N, N), dtype=np.float64)
    fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_BACKWARD', flags=('FFTW_MEASURE',), threads=12)
    y = fft_object(y_hat)
    return y


def func(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    return np.exp(x)*np.sin(y)

dx = 0.02

x = np.arange(-1, 1, dx)
y = np.arange(-1, 1, dx)
X, Y = np.meshgrid(x, y)

N = len(x)

kxw, kyw = np.meshgrid(np.fft.rfftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapw = -4*np.pi**2*(kxw**2+kyw**2)

kxnp, kynp = np.meshgrid(np.fft.fftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapnp = -4*np.pi**2*(kxnp**2+kynp**2)



z = func(X, Y)

lap_z_w = ifftw_(Lapw*fftw_(z))
lap_z_np = np.fft.ifft2(Lapnp*np.fft.fft2(z))

lap_z_np_mag = np.abs(lap_z_np)
lap_z_np_ang = np.angle(lap_z_np)

plt.imshow(z, cmap='plasma')
plt.colorbar()
plt.savefig("f.png", dpi=200)
plt.clf()

plt.imshow(lap_z_w, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_fftw.png", dpi=200)
plt.clf()

plt.imshow(lap_z_np_mag, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_mag.png", dpi=200)
plt.clf()

plt.imshow(lap_z_np_ang, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_ang.png", dpi=200)
plt.clf()

Здесь np.ndarray с именами Lapw и Lapnp - это то, что я думал должен делать дискретный лапласиан. И функция, которую я выбрал, eˣsin (y) , является гармонической функцией, поэтому ее лапласиан должен быть нулевым.

Но результаты программы очень далеки от этого ожидаемого значения. В частности я получаю:

Оригинальная функция f The original funciton

"Лапласиан" f с pyFFTW The

Величина и фаза "лапласиана" f с Numpy The magnitude The phase

Глядя на значения этих графиков (пожалуйста, обратите внимание на диапазон на цветовой панели и тот факт, что 20000 не является каким-либо приличным приближением к 0), становится понятно, почему созданная мной программа дает переполнение, но я не не знаю, как это исправить. Любая помощь будет принята с благодарностью.

1 Ответ

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

Ошибка здесь в том, что ваши предположения о вашей функции неверны.e ^ x sin (y) может показаться гармоничным, но вы рассчитали его только для -1 fft будет периодически неявно продолжать его, т. е. вы получите разрывы на всех границах своей функции.Если функция не является непрерывной, она не является гармонической, и особенно вы получаете расхождения в трансфоме Фурье.Это то, что заставляет ваше БПФ расходиться по краям.Кроме того, «далеко» от краев, результаты выглядят как ожидалось.

...