4

我想使用傅里叶变换在周期性边界条件下找到模拟实体的中心；周期性边界条件意味着，每当某物从盒子的一侧退出时，它就会像经典游​​戏小行星一样扭曲出现在另一侧。

所以我所拥有的是每个时间帧的矩阵（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]) )

会给 在此处输入图像描述

但我不知道如何从这个角度获得峰值的位置。

4

2 回答 2

6

当您只需要一个傅立叶系数时，使用 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

于 2013-08-11T02:34:33.590 回答
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);
于 2014-05-06T12:10:15.603 回答