0

I am trying to compute a 3D power spectrum -- that is, averaged power in frequency shells. I think I'm doing the calculation of shell densities correctly, I'm just not sure how to determine the frequency of each shell.

Supposing the sampling rate Fs is the same in each dimension, and the length of the original samples in each dimension is also the same value N. The shell "index" is idx = sqrt(i*i + j*j + k*k) where i, j, and k are the extent in each direction. How do I compute the frequency of this shell?

cgreen
  • 351
  • 1
  • 17
  • 1
    What is a "frequency shell"? This sounds more like a physics question than a programming question. In general, the discrete Fourier transform of a real-valued spatial-domain volume data set will produce a complex-valued frequency-domain volume data set. It seems like you want something different. – MooseBoys Apr 19 '16 at 20:51
  • Not completely sure what you are trying to ask. Are you trying to ask how to setup an FFTW plan, execute the plan, or interpret the results? – RyanP Apr 19 '16 at 20:53
  • @MooseBoys I mean I want to average the density values (f x f*) for each frequency. It's a shell because it's 3D, it'd be a circle in 2D. – cgreen Apr 19 '16 at 20:53
  • @RyanP sorry, asking for how to interpret the results. After I call `fftw_plan_dft_r2c_3d` I get density values in three dimensions. I average these along each shell (where `radius = sqrt(i*i + j*j + k*k)` is the same). My question is, what is the frequency at that radius? – cgreen Apr 19 '16 at 20:55
  • There are three frequencies, one in each dimension, not a single frequency. – Paul R Apr 19 '16 at 21:04
  • @PaulR true, but the system is symmetric (both in size and sampling rate) and I'm taking the average over a shell. – cgreen Apr 19 '16 at 21:06
  • Well it's easy enough to calculate the frequency in each dimension for a given output point, but I'm not sure that you can talk about a single frequency for one of your spherical "shells". – Paul R Apr 19 '16 at 21:07

1 Answers1

2

It's simpler than you think: the vector (i, j, k) is already the wave vector, so you get the associated frequency by taking its length and dividing it by the length of the edge of your cube.

f = sqrt(i*i + j*j + k*k)/edgelength

The result is the spacial frequency. If you are looking for a temporal frequency, you need some additional information that links the two together.

The only thing that you need to take care of, is the location of the zero frequency within the fft transformed cube: Some algorithms will place it at the upper left corner, others will place it in the center. Wherever it is, you need to take care that you don't misinterpret a low frequency for an aliased high frequency in the opposite direction, i. e. the absolute value of i, j, k should always be less than half the width of your cube.

cmaster - reinstate monica
  • 38,891
  • 9
  • 62
  • 106
  • 1
    Thanks! In [this answer](http://stackoverflow.com/questions/7674877/how-to-get-frequency-from-fft-result) for 1D systems, @Paul R says the frequency of bin `i` is `freq = i * Fs / N`. Do I similarly compute the `unitlength` by dividing the sampling rate by the 1D extent? – cgreen Apr 19 '16 at 21:14
  • 1
    My original version was wrong: the divisor is not the unit length, but the length of the sample, see the updated version of my answer. Paul R's answer is correct: `N/Fs` is the length of the sample, so `i*Fs/N` is the same as `i/samplelength`. – cmaster - reinstate monica Apr 19 '16 at 21:21