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

Question

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

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

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

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

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

• 我对新计算的网格进行 FFT 逆变换。

因此，绘制每个箱上的值的一维直方图将给我我想要的东西。

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

如果我在一个小盒子(300x300x300 cm³ 且点之间间隔至少 1 cm)上运行代码，我不会遇到此问题。不过我确实得到了 600x600x600 cm3 的尺寸。

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

我在这里将我的代码放在一个完整的 block 中，以便轻松复制。

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 次。分解出这样的计算不仅更快，而且产生更可读的代码。