给定一个二维数组，其中每一行代表一个粒子的位置向量，如何有效地计算均方位移（使用fft）？
均方位移定义为
其中r（m）是m行的位置向量，n是行数。
最佳答案
下面的msd方法是直接的，但它是O（N**2）（我改编了这个stackoverflow answer by user morningsun的代码）
def msd_straight_forward(r):
shifts = np.arange(len(r))
msds = np.zeros(shifts.size)
for i, shift in enumerate(shifts):
diffs = r[:-shift if shift else None] - r[shift:]
sqdist = np.square(diffs).sum(axis=1)
msds[i] = sqdist.mean()
return msds
但是，我们可以使用fft使此代码更快。下面的考虑和产生的算法来自this paper，我将只展示如何在python中实现它。
我们可以按以下方式拆分MSD
因此，s_2（m）只是位置的自相关。注意，在一些教科书中，s_（m）表示为自相关（约定a），在一些教科书中，s_（m）*（n-m）表示为自相关（约定b）。
根据Wiener-Khinchin定理，函数的功率谱密度（PSD）是自相关的Fourier变换。
这意味着我们可以计算一个信号的psd，然后对其进行傅里叶逆变换，得到自相关（在惯例b中）。对于离散信号，我们得到循环自相关。
然而，通过对数据进行零填充，我们可以得到非循环自相关算法如下
def autocorrFFT(x):
N=len(x)
F = np.fft.fft(x, n=2*N) #2*N because of zero-padding
PSD = F * F.conjugate()
res = np.fft.ifft(PSD)
res= (res[:N]).real #now we have the autocorrelation in convention B
n=N*np.ones(N)-np.arange(0,N) #divide res(m) by (N-m)
return res/n #this is the autocorrelation in convention A
对于术语S 1（m），我们利用这样一个事实，即可以找到（N-m）*s1（m）的递归关系（这在第4.2节中的paper中进行了解释）。
我们定义
并通过
这将产生以下均方位移的代码
def msd_fft(r):
N=len(r)
D=np.square(r).sum(axis=1)
D=np.append(D,0)
S2=sum([autocorrFFT(r[:, i]) for i in range(r.shape[1])])
Q=2*D.sum()
S1=np.zeros(N)
for m in range(N):
Q=Q-D[m-1]-D[N-m]
S1[m]=Q/(N-m)
return S1-2*S2
您可以比较msd_straight_forward（）和msd_fft（），会发现它们产生相同的结果，尽管msd_fft（）对于大型N
一个小基准：用
r = np.cumsum(np.random.choice([-1., 0., 1.], size=(N, 3)), axis=0)
对于n=100000，我们得到
$ %timeit msd_straight_forward(r)
1 loops, best of 3: 2min 1s per loop
$ %timeit msd_fft(r)
10 loops, best of 3: 253 ms per loop