时间戳

首页 > 解决方案 > 傅里叶空间中的滤波器的行为不像它应该的那样

python - 傅里叶空间中的滤波器的行为不像它应该的那样

问题描述

这是我提出的已回答问题的后续行动，可在此处找到。

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

根据我在傅里叶空间中使用滤波器的解析表达式遇到的问题，我现在在真实空间中生成球并采用它的 FFT 来获得我的滤波器。但是，我从这种方法中得到的直方图真的很奇怪，我希望得到一个高斯分布： 我的盒子中质量密度的归一化直方图

我的代码如下：

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 中所占的体积减小，该值减小。我检查了它是否给了我想要的半径的球体。可以在此处找到有关该主题的详细信息（第 2.5 部分和图 8 以确保准确性）。

- 编辑 -

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

标签： pythonnumpyfilteringfft

解决方案

我的问题来自我如何生成过滤器。在我的代码中，我将权重与不完全在球体中的体素相关联的方式是不连续的：例如，我将权重 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

推荐阅读