2

I am trying to calculate the power spectrum (1D) of a regular distributed 3D density database. But I keep getting the peak of power spectrum at 0. I am appreciating for your help!

Here is my code:

import numpy as np
import matplotlib.pyplot as plt

bins = 1.
tbins = 400.

density = np.random.rand(400,400,400)

x, y, z = np.mgrid[0:400, 0:400, 0:400]
x = x - 199.5
y = y - 199.5
z = z - 199.5

dist = np.sqrt(x**2+y**2+z**2)

FT = np.fft.fftn(density)
power = FT.real*FT.real + FT.imag*FT.imag

P = power.reshape(np.size(power))
dist = dist.reshape(np.size(dist))

intervals = np.array([nn*bins for nn in range(0,int(tbins)+1)])

p = np.histogram(dist, bins=intervals, weights=P)[0]
pd = np.histogram(dist, bins=intervals)[0]
pd.astype('float')
p = p/pd

plt.figure()
plt.plot(2.*np.pi/intervals[1:], p)
plt.show()
rtime
  • 349
  • 6
  • 23
Donna
  • 21
  • 2
  • 1
    If your data has a non-zero mean then you would expect a large peak at 0. – Paul R Feb 10 '17 at 03:37
  • 1
    You’re using `rand`, i.e., uniform between 0 and 1, which means `density` has a large non-zero mean. So as @PaulR mentioned, you’ll get an enormous `FT[0,0,0]` compared to every other element of `FT`. Try `randn` (note the `n` at the end)—the peak at 0-bin will be much smaller since `density` will be zero-mean Gaussian. – Ahmed Fasih Feb 10 '17 at 04:10

0 Answers0