使用python和FFT计算均方位移

给定一个二维数组,其中每一行代表一个粒子的位置矢量,如何有效地计算均方位移(使用FFT)？均方位移定义为

其中r(m)是行m的位置矢量,N是行数.

以下为msd的直接方法有效,但它是O(N**2)(我使用用户morningsun修改了此stackoverflow应答的代码)

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更快地编写代码.下面的考虑和由此产生的算法来自本文,我将展示如何在python中实现它.我们可以通过以下方式拆分MSD

因此,S_2(m)只是位置的自相关.注意,在一些教科书中,S_2(m)表示为自相关(约定A),而在一些S_2(m)*(Nm)中表示为自相关(约定B).根据Wiener-Khinchin定理,函数的功率谱密度(PSD)是自相关的傅里​​叶变换.这意味着我们可以计算信号的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),我们利用了这样一个事实:可以找到(Nm)*S_1(m)的递归关系(本文在4.2节中对此进行了解释).我们定义

并找到S_1(m)via

这产生以下均方位移代码

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 = 100.000,我们得到

$ %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