title: Power Spectrum and Autocorrelation of Data in Numpy
body: <p>I am interested in computing the power spectrum of a system of particles (~100,000) in 3D space with Python. What I have found so far is a group of functions in Numpy (<code>fft</code>,<code>fftn</code>,..) which compute the discrete Fourier transform, of which the square of the absolute value is the power spectrum. My question is a matter of how my data are being represented - and truthfully may be fairly simple to answer. </p>
  
  <p>The data structure I have is an array which has a shape of (<em>n</em>,2), <em>n</em> being the number of particles I have, and each column representing either the x, y, and z coordinate of the <em>n</em> particles. The function I believe I should be using it the <code>fftn()</code> function, which takes the discrete Fourier transform of an <em>n</em>-dimensional array - but it says nothing about the format. How should the data be represented as a data structure to be fed into <code>fftn</code>?</p>
  
  <p>Here is what I've tried so far to test the function:</p>
  
  <pre><code>import numpy as np
  import random
  import matplotlib.pyplot as plt
  
  DATA = np.zeros((100,3))
  
  for i in range(len(DATA)):
      DATA[i,0] = random.uniform(-1,1)
      DATA[i,1] = random.uniform(-1,1)
      DATA[i,2] = random.uniform(-1,1)
  
  FFT = np.fft.fftn(DATA)
  PS = abs(FFT)**2
  
  plt.plot(PS)
  plt.show()
  </code></pre>
  
  <p>The array entitled <code>DATA</code> is a mock array, ultimately the thing which will be 100,000 by 3 in shape. The output of the code gives me something like:
  <img src="https://i.stack.imgur.com/k0ey3.png" alt="enter image description here"></p>
  
  <p>As you can see, I think this is giving me three 1D power spectra (1 for each column of my data), but really I'd like a power spectrum as a function of radius.</p>
  
  <p>Does anybody have any advice or alternative methods/packages they know of to compute the power spectrum (I'd even settle for the two point autocorrelation function).</p>