Using FFT to find the center of mass under periodic boundary conditions

Question

I would like to use the Fourier transform to find the center of a simulated entity under periodic boundary condition; periodic boundary conditions means, that whenever something exits through one side of the box, it is warped around to appear on the opposite side just like in the classic game asteroids.

So what I have is for each time frame a matrix (Nx3) with N the number of points in xyz. what I want to do is determine the center of that cloud even if it all moved over the periodic boundary and is so to say stuck in between.

My idea for an solution would now be do a (mass weigted) histogram of these points and then perform an FFT on that and use the phase of the first Fourier coefficient to determine where in the box the maximum would be.

as a test case I have used

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] ]

So now I can use FFT since it expects a periodic function anyways.

## 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)

That way I get a sine wave with its maximum exactly where the maximum of the histogram is.

Now I'd like to extract the position of the maximum not by taking the max function on the sine vector, but somehow it should be retrievable from the first (not the 0th) Fourier coefficient, since that should somehow contain the phase shift of the sine to have its maximum exactly at the maximum of the histogram.

Indeed, plotting

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

will give

But I can't figure out how to get the position of the peak from this angle.

Answer 1

Using the FFT is overkill when all you need is one Fourier coefficent. Instead, you can simply compute the dot product of your data with

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

where N is the number of points. (The time to compute all the Fourier coefficients with the FFT is O(N*log(N)). Computing just one coefficient is O(N).)

Here's a script similar to yours. The data is put in y; the coordinates of the data points are in 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()

This generates the plot:

To answer the question "Why F1.conj()?":

The complex conjugate of F1 is used because of the minus sign in w = np.exp(-2j*np.pi*np.arange(N) / N) (which I used because it is a common convention).

Since w can be written

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)

the dot product y.dot(w) is basically a projection of y onto cos(2*pi*arange(N)/N) (the real part of F1) and -sin(2*pi*arange(N)/N) (the imaginary part of F1). But when we figure out the phase of the maximum, it is based on the functions cos(...) and sin(...). Taking the complex conjugate accounts for the opposite sign of the sin() function. If w = np.exp(2j*np.pi*np.arange(N) / N) were used instead, the complex conjugate of F1 would not be needed.

Answer 2

You could calculate the circular mean directly on your data.

When calculating the circular mean, your data is mapped to -pi..pi. This mapped data is interpreted as angle to a point on the unit circle. Then the mean value of x and y component is calculated. The next step is to calculate the resulting angle and map it back to the defined "box".

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);