julia - FFTW.jl 用于 2D 阵列：扩散只发生在 1D

Question

根据我的阅读，当 A 是 2D 数组时，使用 FFTW.jl / AbstractFFTs.jl 的 fft(A) 应该在 2D 中执行 fft，而不是按列执行。知道为什么当（我认为）我将缩放的二阶空间导数添加到 u(t,x) 时，我只看到列方向的扩散，就好像使用显式求解器一样？

谢谢！我对此很陌生。

代码和热图截图

using Random
using FFTW
using Plots

gr()

N = (100,100)

# initialize with gaussian noise
u = randn(Float16, (N[1], N[2])).*0.4.+0.4
# include square of high concentration to observe diffusion clearly
u[40:50,40:50] .= 3

N = size(x)
L = 100
k1 = fftfreq(51)
k2 = fftfreq(51)
lap_mat = -(k1.^2 + k2.^2)

function lap_fft(x)
    lapF = rfft(x)
    lap = irfft(lap_mat.*lapF, 100)
    return lap
end

# ode stepper or Implicit-Explicit solver
for i in 1:100000
    u+=lap_fft(u)*0.0001
end
# plot state
heatmap(u)

Answer 1

仅仅因为您正在执行真正的 FFT，并不意味着您可以对结果进行真正的逆 fft。rfft 来自 R -> C。但是，您可以执行以下操作：

function lap_fft(x)
    lapF = complex(zeros(100,100));        # only upper half filled
    lapF[1:51,1:100] = rfft(x) .* lap_mat; # R -> C
    return abs.(ifft(lapF));               # C -> R
end

真正的 FFT 到复杂的频域（由于数据冗余，只填充了上半部分），在频域中乘以你的滤波器，将 FFT 逆向复杂的图像域并获得幅度abs.()，实部real.()等。
但老实说，为什么要使用真正的 fft 呢？

using Random
using FFTW
using Plots

gr()

N = (100,100)

# initialize with gaussian noise
u = randn(Float16, (N[1], N[2])).*0.4.+0.4;
# include square of high concentration to observe diffusion clearly
u[40:50,40:50] .= 3;

N = size(u);
L = 100;
k1 = fftfreq(100);
k2 = fftfreq(100);
tmp = -(k1.^2 + k2.^2);
lap_mat = sqrt.(tmp.*reshape(tmp,1,100));


function lap_fft(x)
    return abs.(ifftshift(ifft(fftshift(ifftshift(fft(fftshift(x))).*lap_mat))));
end

# ode stepper or Implicit-Explicit solver
for i in 1:100000
    u+=lap_fft(u)*0.001;
end
# plot state
heatmap(u)