逆FFT不应该返回负值 [英] Inverse FFT returns negative values when it should not
问题描述
我在3D盒子中有几个点(x，y，z坐标)，并带有相关的质量.我想绘制在给定半径R的球体中发现的质量密度的直方图.
I have several points (x,y,z coordinates) in a 3D box with associated masses. I want to draw an histogram of the mass-density that is found in spheres of a given radius R.
如果没有出现我认为可能有的任何错误，我编写的代码将以下列方式工作:
I have written a code that, providing I did not make any errors which I think I may have, works in the following way:
-
我的真实"数据是巨大的，因此我写了一些代码来随机产生任意重叠的非重叠点.
My "real" data is something huge thus I wrote a little code to generate non overlapping points randomly with arbitrary mass in a box.
我计算出3D直方图(按质量加权)，其装仓比我的球体半径小10倍.
I compute a 3D histogram (weighted by mass) with a binning about 10 times smaller than the radius of my spheres.
我对直方图进行FFT，计算波模(kx，ky和kz)，并使用它们通过3D top的解析表达式将傅立叶空间中的直方图相乘傅里叶空间中的-hat窗口(球面过滤)功能.
I take the FFT of my histogram, compute the wave-modes (kx, ky and kz) and use them to multiply my histogram in Fourier space by the analytic expression of the 3D top-hat window (sphere filtering) function in Fourier space.
我对新计算的网格进行了逆FFT.
I inverse FFT my newly computed grid.
因此在每个bin上绘制值的一维直方图会给我我想要的东西.
Thus drawing a 1D-histogram of the values on each bin would give me what I want.
我的问题如下:鉴于我所做的事情，我的反向FFT网格中应该没有任何负值(第4步)，但是我得到了一些负值，并且其值比数值误差高得多.
My issue is the following: given what I do there should not be any negative values in my inverted FFT grid (step 4), but I get some, and with values much higher that the numerical error.
如果我在一个小盒子(300x300x300 cm 3 )上运行我的代码，并且这些点之间的距离至少为1 cm，则不会出现此问题.我确实以600x600x600 cm3的尺寸得到它.
If I run my code on a small box (300x300x300 cm3 and the points of separated by at least 1 cm) I do not get the issue. I do get it for 600x600x600 cm3 though.
如果将所有质量都设置为0，从而在一个空网格上工作，那么我确实会得到0，而没有任何明显的问题.
If I set all the masses to 0, thus working on an empty grid, I do get back my 0 without any noted issues.
我在这里以完整的块形式提供我的代码，以便将其轻松复制.
I here give my code in a full block so that it is easily copied.
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()
您知道为什么我会得到这些负值吗?您是否认为有一种更简单的方法?
Do you know why I'm getting these negative values? Do you think there is a simpler way to proceed?
很抱歉，如果这是一篇很长的文章，我试图让它变得很清楚，并会根据您的评论进行编辑，非常感谢！
Sorry if this is a very long post, I tried to make it very clear and will edit it with your comments, thanks a lot!
-编辑-
A follow-up question on the issue can be found [here].1
推荐答案
您在频域中创建的滤波器只是您要创建的滤波器的近似值.问题在于我们在这里处理DFT，而不是连续域FT(具有无限频率).球的傅立叶变换确实是您所描述的函数，但是此函数无限大-它不受频带限制！
The filter you create in the frequency domain is only an approximation to the filter you want to create. The problem is that we are dealing with the DFT here, not the continuous-domain FT (with its infinite frequencies). The Fourier transform of a ball is indeed the function you describe, however this function is infinitely large -- it is not band-limited!
通过仅在窗口中对该函数进行采样，可以有效地将其与理想的低通滤波器(域的矩形)相乘.在空间域中，该低通滤波器具有负值.因此，您创建的过滤器在空间域中也具有负值.
By sampling this function only within a window, you are effectively multiplying it with an ideal low-pass filter (the rectangle of the domain). This low-pass filter, in the spatial domain, has negative values. Therefore, the filter you create also has negative values in the spatial domain.
这是通过Kh的逆变换的原点的切片(在我应用fftshift之后将原点移动到图像的中间，以便更好地显示):
This is a slice through the origin of the inverse transform of Kh (after I applied fftshift to move the origin to the middle of the image, for better display):
您可以在这里说到，有些振铃会导致负值.
As you can tell here, there is some ringing that leads to negative values.
克服振铃的一种方法是在频域中应用开窗功能.另一种选择是在空间域中生成一个球，然后计算其傅立叶变换.第二种选择将是最简单的实现.请记住，在空间域中的内核还必须在左上像素处具有原点，以获取正确的FFT.
One way to overcome this ringing is to apply a windowing function in the frequency domain. Another option is to generate a ball in the spatial domain, and compute its Fourier transform. This second option would be the simplest to achieve. Do remember that the kernel in the spatial domain must also have the origin at the top-left pixel to obtain a correct FFT.
通常在空间域中应用开窗函数，以避免在计算FFT时出现图像边界问题.在这里，我建议在频域中应用这样的窗口，以避免在计算IFFT时出现类似的问题.但是请注意，这将始终进一步减少内核的带宽(毕竟，开窗功能将充当低通滤波器)，因此在空间域(即空间域)中产生前景到背景的平滑过渡内核将没有您想要的那样尖锐的过渡).最著名的窗口功能是 Hamming和Hann窗口，但是有很多其他人值得一试.
A windowing function is typically applied in the spatial domain to avoid issues with the image border when computing the FFT. Here, I propose to apply such a window in the frequency domain to avoid similar issues when computing the IFFT. Note, however, that this will always further reduce the bandwidth of the kernel (the windowing function would work as a low-pass filter after all), and therefore yield a smoother transition of foreground to background in the spatial domain (i.e. the spatial domain kernel will not have as sharp a transition as you might like). The best known windowing functions are Hamming and Hann windows, but there are many others worth trying out.
我简化了将Kh计算为以下代码的代码:
I simplified your code to compute Kh to the following:
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次.排除这种计算不仅速度更快，而且产生了更具可读性的代码.
I find this easier to read than the nested loops. It should also be significantly faster, and avoid the need for njit. Note that you were computing the same distance (what I call kr here) 5 times. Factoring out such computation is not only faster, but yields more readable code.
这篇关于逆FFT不应该返回负值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
