使用FFT在周期性边界条件下找到质心

Question

我想用傅立叶变换在周期边界条件下找到模拟实体的中心; 周期性边界条件意味着,只要有东西从盒子的一侧出来,就会像在经典游戏小行星中一样扭曲出现在对面.

所以我所拥有的是每个时间帧的矩阵(Nx3),其中N是xyz中的点数.我想要做的是确定云的中心,即使它们都在周期性边界上移动,也就是说卡在两者之间.

我对解决方案的想法现在是对这些点进行(质量权重)直方图,然后对其进行FFT并使用第一个傅里叶系数的相位来确定最大值在框中的位置.

作为我用过的测试用例

import numpy as np
Points_x  = np.random.randn(10000)
Box_min   = -10
Box_max   =  10
X         = np.linspace( Box_min, Box_max, 100 )

### make a Histogram of the points
Histogram_Points = np.bincount( np.digitize( Points_x, X ),  minlength=100 )

### make an artifical shift over the periodic boundary
Histogram_Points = np.r_[ Histogram_Points[45:], Histogram_Points[:45] ]

所以现在我可以使用FFT,因为它总是需要一个周期函数.

## doing fft
F = np.fft.fft(Histogram_Points)

## getting rid of everything but first harmonic
F[2:] = 0.

## back transforming
Fist_harmonic = np.fft.ifft(F)

这样我得到的正弦波的最大值恰好是直方图的最大值.

现在我想提取最大值的位置不是通过在正弦向量上取最大函数,但不知何故它应该可以从第一个(不是第0个)傅立叶系数中检索,因为它应该以某种方式包含相移正弦最大值恰好在直方图的最大值.

的确,密谋

Cos_approx = cos( linspace(0,2*pi,100) * angle(F[1]) )

会给

但我无法弄清楚如何从这个角度获得峰值的位置.

Answer 1

当你需要的只是一个傅里叶系数时,使用FFT是过度的.相反,您可以简单地计算数据的点积

w = np.exp(-2j*np.pi*np.arange(N) / N)

其中N是点数.(用FFT计算所有傅立叶系数的时间是O(N*log(N)).只计算一个系数是O(N).)

这是一个类似于你的脚本.数据被放入y; 数据点的坐标在x.

import numpy as np

N = 100

# x coordinates of the data
xmin = -10
xmax = 10
x = np.linspace(xmin, xmax, N, endpoint=False)

# Generate data in y.
n = 35
y = np.zeros(N)
y[:n] = 1 - np.cos(np.linspace(0, 2*np.pi, n))
y[:n] /= 0.7 + 0.3*np.random.rand(n)
m = 10
y = np.r_[y[m:], y[:m]]

# Compute coefficent 1 of the discrete Fourier transform.
w = np.exp(-2j*np.pi*np.arange(N) / N)
F1 = y.dot(w)
print "F1 =", F1

# Get the angle of F1 (in the interval [0,2*pi]).
angle = np.angle(F1.conj())
if angle < 0:
    angle += 2*np.pi

center_x = xmin + (xmax - xmin) * angle / (2*np.pi)
print "center_x = ", center_x

# Create the first sinusoidal mode for the plot.
mode1 = (F1.real * np.cos(2*np.pi*np.arange(N)/N) -
         F1.imag*np.sin(2*np.pi*np.arange(N)/N))/np.abs(F1)


import matplotlib.pyplot as plt

plt.clf()
plt.plot(x, y)
plt.plot(x, mode1)
plt.axvline(center_x, color='r', linewidth=1)
plt.show()

这会生成情节:

要回答"为什么F1.conj()？" 这个问题:

F1由于减号w = np.exp(-2j*np.pi*np.arange(N) / N)(我使用因为它是一种常见的惯例)而使用 复共轭.

既然w可以写

w = np.exp(-2j*np.pi*np.arange(N) / N)
  = cos(-2*pi*arange(N)/N) + 1j*sin(-2*pi*arange(N)/N)
  = cos(2*pi*arange(N)/N) - 1j*sin(2*pi*arange(N)/N)

点积y.dot(w)基本上是一个投影y到 cos(2*pi*arange(N)/N)(的实部F1)和-sin(2*pi*arange(N)/N) (的虚部F1).但是当我们计算出最大值的相位时,它基于函数cos(...)和sin(...).采用复共轭解释了sin()函数的相反符号.如果w = np.exp(2j*np.pi*np.arange(N) / N)使用,则F1不需要复共轭.