傅立叶空间中的滤波器的行为与预期不同

这是我提出的已回答问题的后续内容，可以找到here https://stackoverflow.com/questions/54022376/inverse-fft-returns-negative-values-when-it-should-not.

我在 3D 框中有几个具有关联质量的点（x、y、z 坐标）。我想绘制给定半径球体中质量密度的直方图R。这个想法是计算我的盒子的 3D 直方图（分箱比半径小得多），对其进行 FFT，乘以滤波器（真实空间中的球）并对结果进行逆 FFT。从那里，我只计算每个 3D-bin 中获得的值的 1D 直方图。

Following the issue I had by using an analytic expression of the filter in Fourier space, I am now generating the ball in real space and take its FFT to obtain my filter. However, the histogram I get out of this method is really strange, where I would expect a Gaussian I am getting this:

我的代码如下：

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

size = 100

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


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(1)
    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
R = 10

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

Fh = np.fft.fftn(H, axes=(-3,-2, -1))

# Generate the filter in real space

Kreal = np.zeros((b,b,b))


X = edges[0]
Y = edges[1]
Z = edges[2]

mid = int(b/2)

s = (X.max()-X.min()+Y.max()-Y.min()+Z.max()-Z.min())/(3*b)

cst = 1/2 + (1/12 - (R/s)**2)*np.arctan((0.5*np.sqrt((R/s)**2-0.5))/(0.5-(R/s)**2)) + 1/3*np.sqrt((R/s)**2-0.5) + ((R/s)**2 - 1/12)*np.arctan(0.5/(np.sqrt((R/s)**2-0.5))) - 4/3*(R/s)**3*np.arctan(0.25/((R/s)*np.sqrt((R/s)**2-0.5)))


@njit(parallel=True)
def remp(Kreal):
    for i in range(b):
        for j in range(b):
            for k in range(b):
                a = cst - np.sqrt((X[i]-X[mid])**2 + (Y[j]-Y[mid])**2 + (Z[k]-Z[mid])**2)/s
                if a >= 0.1 and a < 0.2:
                    Kreal[i][j][k] = 0.1
                elif a >= 0.2 and a < 0.3:
                   Kreal[i][j][k] = 0.2 
                elif a >= 0.3 and a < 0.4:
                   Kreal[i][j][k] = 0.3
                elif a >= 0.4 and a < 0.5:
                   Kreal[i][j][k] = 0.4
                elif a >= 0.5 and a < 0.6:
                   Kreal[i][j][k] = 0.5
                elif a >= 0.6 and a < 0.7:
                   Kreal[i][j][k] = 0.6
                elif a >= 0.7 and a < 0.8:
                   Kreal[i][j][k] = 0.7
                elif a >= 0.8 and a < 0.9:
                   Kreal[i][j][k] = 0.8
                elif a >= 0.9 and a < 0.99:
                   Kreal[i][j][k] = 0.9
                elif a >= 0.99:
                   Kreal[i][j][k] = 1
    return Kreal


Kreal = remp(Kreal)

Kreal = np.fft.ifftshift(Kreal)

Kh = np.fft.fftn(Kreal, axes=(-3,-2, -1))

Gh = np.multiply(Fh, Kh)

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

# Generate the filter in fourier space using its analytic expression

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))

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


Gh = np.multiply(Fh, Kh)

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


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

D2 = Density2.flatten()
N2 = np.mean(D2)


# 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,'.',label = "Defining the Filter in real space")


hist, bins = np.histogram(D2/N2, bins='auto', density=True)
bin_centers = (bins[1:]+bins[:-1])*0.5
plt.plot(bin_centers, hist,'.',label = "Using analytic expression")



plt.xlabel('Normalised Density')
plt.ylabel('Probability density')
plt.legend()
plt.show()

你明白为什么会发生这种情况吗？非常感谢您的帮助。

PS：一长串if当我在真实空间中定义过滤器时的语句来自我在网格上绘制球体的方式。我将值 1 分配给所有 100% 在球体内的 bin，然后该值随着 bin 中球体占据的体积减小而减小。我检查它给了我一个所需半径的球体。有关该主题的详细信息可以找到here http://geonumerics.mit.edu/assets/docs/papers/dem_lbm_mapping.pdf（第 2.5 部分和图 8 的准确性）。

--EDIT--

仅当所有粒子质量相同时，代码才会表现得像这样

我的问题来自于我如何生成过滤器。在我的代码中，我将权重与不完全位于球体中的体素相关联的方式是不连续的：例如，我将权重 0.1 赋予体积比在 0.1 和 0.2 之间的体素。

因此，当所有点具有相同的质量时，会发生什么：我的网格中有 1 的倍数，我将其与有限数量的系数相乘，因此我的网格可以采用有限数量的可能值，因此有些箱是空的或至少“不太满”。当粒子的质量分布更加连续时，这种情况就不太可能发生。

因此，解决方法是为体素指定正确的权重。

def remp(Kreal):
for i in range(b):
    for j in range(b):
        for k in range(b):
            a = cst - np.sqrt((X[i]-X[mid])**2 + (Y[j]-Y[mid])**2 + (Z[k]-Z[mid])**2)/s
            if a >= 0.1 and a < 0.99:
                Kreal[i][j][k] = a
            elif a >= 0.99:
               Kreal[i][j][k] = 1