python - 逆向 FFT 在不应返回负值时返回负值

Question

我在具有相关质量的 3D 框中有几个点(x、y、z 坐标)。我想绘制在给定半径 R 的球体中发现的质量密度直方图。

我已经编写了一个代码，如果我没有犯任何我认为我可能有的错误，它的工作方式如下:

• 我的“真实”数据非常庞大，因此我编写了一些代码来随机生成盒子中任意质量的非重叠点。

• 我计算了一个 3D 直方图(按质量加权)，分箱比我的球体半径小 10 倍。

• 我对直方图进行 FFT，计算波形模式(kx、ky 和 kz)并使用它们将傅立叶空间中的直方图乘以傅立叶空间中的 3D 顶帽窗口(球体滤波)函数的解析表达式。

• 我对新计算的网格进行 FFT 反演。

因此绘制每个 bin 值的一维直方图会给我想要的结果。

我的问题如下:考虑到我所做的事情，我的反向 FFT 网格(第 4 步)中不应该有任何负值，但我得到了一些，并且值比数值误差高得多。

如果我在一个小盒子上运行我的代码(300x300x300 cm³ 并且点之间至少相隔 1 cm)我没有得到这个问题。不过我确实得到了 600x600x600 立方厘米。

如果我将所有质量设置为 0，从而在空网格上工作，我确实会返回我的 0 而不会出现任何值得注意的问题。

我在这里完整地给出了我的代码，以便于复制。

import numpy as np
import matplotlib.pyplot as plt
import random
from numba import njit

# 1. Generate a bunch of points with masses from 1 to 3 separated by a radius of 1 cm

radius = 1
rangeX = (0, 100)
rangeY = (0, 100)
rangeZ = (0, 100)
rangem = (1,3)
qty = 20000  # or however many points you want

# Generate a set of all points within 1 of the origin, to be used as offsets later
deltas = set()
for x in range(-radius, radius+1):
    for y in range(-radius, radius+1):
        for z in range(-radius, radius+1):
            if x*x + y*y + z*z<= radius*radius:
                deltas.add((x,y,z))

X = []
Y = []
Z = []
M = []
excluded = set()
for i in range(qty):
    x = random.randrange(*rangeX)
    y = random.randrange(*rangeY)
    z = random.randrange(*rangeZ)
    m = random.uniform(*rangem)
    if (x,y,z) in excluded: continue
    X.append(x)
    Y.append(y)
    Z.append(z)
    M.append(m)
    excluded.update((x+dx, y+dy, z+dz) for (dx,dy,dz) in deltas)

print("There is ",len(X)," points in the box")

# Compute the 3D histogram
a = np.vstack((X, Y, Z)).T
b = 200

H, edges = np.histogramdd(a, weights=M, bins = b)      

# Compute the FFT of the grid
Fh = np.fft.fftn(H, axes=(-3,-2, -1))

# Compute the different wave-modes
kx = 2*np.pi*np.fft.fftfreq(len(edges[0][:-1]))*len(edges[0][:-1])/(np.amax(X)-np.amin(X))
ky = 2*np.pi*np.fft.fftfreq(len(edges[1][:-1]))*len(edges[1][:-1])/(np.amax(Y)-np.amin(Y))
kz = 2*np.pi*np.fft.fftfreq(len(edges[2][:-1]))*len(edges[2][:-1])/(np.amax(Z)-np.amin(Z))

# I create a matrix containing the values of the filter in each point of the grid in Fourier space

R = 5                                                                                               
Kh = np.empty((len(kx),len(ky),len(kz)))

@njit(parallel=True)
def func_njit(kx, ky, kz, Kh):
    for i in range(len(kx)):
        for j in range(len(ky)):
            for k in range(len(kz)):
                if np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2) != 0:
                    Kh[i][j][k] = (np.sin((np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2))*R)-(np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2))*R*np.cos((np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2))*R))*3/((np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2))*R)**3
                else:
                    Kh[i][j][k] = 1
    return Kh

Kh = func_njit(kx, ky, kz, Kh)

# I multiply each point of my grid by the associated value of the filter (multiplication in Fourier space = convolution in real space)

Gh = np.multiply(Fh, Kh)

# I take the inverse FFT of my filtered grid. I take the real part to get back floats but there should only be zeros for the imaginary part.

Density = np.real(np.fft.ifftn(Gh,axes=(-3,-2, -1)))

# Here it shows if there are negative values the magnitude of the error

print(np.min(Density))

D = Density.flatten()
N = np.mean(D)

# I then compute the histogram I want

hist, bins = np.histogram(D/N, bins='auto', density=True)
bin_centers = (bins[1:]+bins[:-1])*0.5
plt.plot(bin_centers, hist)
plt.xlabel('rho/rhom')
plt.ylabel('P(rho)')

plt.show()

你知道我为什么会得到这些负值吗？您认为有更简单的方法可以继续吗？

很抱歉，如果这是一篇很长的帖子，我已经尽量把它讲清楚了，并且会根据您的意见对其进行编辑，非常感谢!

-编辑-

可以在 [此处] 找到有关该问题的后续问题。 1

Answer 1

您在频域中创建的滤波器只是您要创建的滤波器的近似值。问题是我们在这里处理的是 DFT，而不是连续域 FT(具有无限频率)。球的傅里叶变换确实是你描述的函数，但是这个函数是无限大的——它不是带限的!

通过仅在窗口内对该函数进行采样，您可以有效地将其与理想的低通滤波器(域的矩形)相乘。该低通滤波器在空间域中具有负值。因此，您创建的过滤器在空间域中也具有负值。

这是 Kh 逆变换原点的切片(在我应用 fftshift 将原点移动到图像中间，以便更好地显示之后) :

正如您在这里所看到的，有一些振铃会导致负值。

克服这种振铃的一种方法是在频域中应用窗口函数。另一种选择是在空间域中生成球，并计算其傅里叶变换。第二种选择是最容易实现的。请记住，空间域中的内核还必须以左上角像素为原点才能获得正确的 FFT。

窗口函数通常应用于空间域，以避免在计算 FFT 时出现图像边界问题。在这里，我建议在频域中应用这样一个窗口，以避免在计算 IFFT 时出现类似问题。但是请注意，这将始终进一步降低内核的带宽(毕竟窗口函数将用作低通滤波器)，因此会在空间域(即空间域)中产生前景到背景的更平滑过渡内核不会像您希望的那样急剧过渡)。最著名的窗口函数是 Hamming and Hann windows , 但有 many others值得一试。

---

不请自来的建议:

我将计算 Kh 的代码简化为以下内容:

kr = np.sqrt(kx[:,None,None]**2 + ky[None,:,None]**2 + kz[None,None,:]**2)
kr *= R
Kh = (np.sin(kr)-kr*np.cos(kr))*3/(kr)**3
Kh[0,0,0] = 1

我发现这比嵌套循环更容易阅读。它也应该明显更快，并且避免需要 njit。请注意，您正在计算相同的距离(我在这里称之为 kr)5 次。分解出这样的计算不仅速度更快，而且会产生更易读的代码。