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

这样我得到一个正弦波，其最大值正好在直方图的最大值处。

现在我想提取最大值的位置，而不是通过在正弦向量上取 max 函数，但不知何故，它应该可以从第一个（而不是第 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 的复共轭。

Answer 2

您可以直接根据数据计算循环平均值。

在计算循环平均值时，您的数据将映射到 -pi..pi。该映射数据被解释为与单位圆上某个点的角度。然后计算 x 和 y 分量的平均值。下一步是计算结果角度并将其映射回定义的“框”。

import numpy as np
import matplotlib.pyplot as plt

Points_x  = np.random.randn(10000)+1
Box_min   = -10
Box_max   =  10
Box_width = Box_max - Box_min

#Maps Points to Box_min ... Box_max with periodic boundaries
Points_x = (Points_x%Box_width + Box_min)
#Map Points to -pi..pi
Points_map = (Points_x - Box_min)/Box_width*2*np.pi-np.pi
#Calc circular mean
Pmean_map  = np.arctan2(np.sin(Points_map).mean() , np.cos(Points_map).mean())
#Map back
Pmean = (Pmean_map+np.pi)/(2*np.pi) * Box_width + Box_min

#Plotting the result
plt.figure(figsize=(10,3))
plt.subplot(121)
plt.hist(Points_x, 100);
plt.plot([Pmean, Pmean], [0, 1000], c='r', lw=3, alpha=0.5);
plt.subplot(122,aspect='equal')
plt.plot(np.cos(Points_map), np.sin(Points_map), '.');
plt.ylim([-1, 1])
plt.xlim([-1, 1])
plt.grid()
plt.plot([0, np.cos(Pmean_map)], [0, np.sin(Pmean_map)], c='r', lw=3, alpha=0.5);