逆 FFT 在不应该返回负值时返回

Question

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

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

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

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

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

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

因此，在每个 bin 上绘制一个一维直方图会给我想要的。

我的问题如下：考虑到我所做的事情，我的倒置 FFT 网格（第 4 步）中不应该有任何负值，但我得到了一些，并且值远高于数值误差。

如果我在一个小盒子上运行我的代码（300x300x300 厘米³ 并且点至少相隔 1 厘米）我不明白这个问题。不过，我确实得到了 600x600x600 cm3。

如果我将所有质量设置为 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()

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

对不起，如果这是一个很长的帖子，我试图把它说得很清楚，并会与你的 cmets 一起编辑它，非常感谢！

-编辑-

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

Answer 1

只是猜测：

您从哪里得到虚部必须为零的想法？你有没有试过取绝对值 (sqrt(re^2 + im^2)) 而忘记相位而不是只取实数？只是想到了一些事情。

Answer 2

您在频域中创建的滤波器只是您要创建的滤波器的近似值。问题是我们在这里处理的是 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 次。分解出这样的计算不仅更快，而且产生更易读的代码。