Forums

FFTW Power Spectrum

Started by Mosby February 5, 2010
Dear all,

for quite a while I struggle with the following problem : 
Having a 3-D Array of equally data in X,Y & Z , I want to calculate the
power spectra. 

=============

#include "fftw3.h"
#include <math.h>
#include <iostream>
#include <complex>
#include <blitz/array.h>

using namespace blitz;

main(int argc,char **argv)
{

 // Blitz++ is an c++ class for faciliating array handling
blitz::Array<double, 3> phi(32,32,32);        // creates 3D Array with
indices 0-31 for x, y, z
blitz::Array<std::complex<double>, 3> phi_k(32,32,17);

fftw_plan plan_Phi3DForward      = fftw_plan_dft_r2c_3d(32, 32, 32, 
phi.data(), (fftw_complex *)  phi_k.data(), FFTW_ESTIMATE);

// 3d symmetric  sin wave
 for(int i = 0; i <= 31; i++) {
 for(int j = 0; j <= 31; j++) {
  for(int k = 0; k <= 31; k++) {
        phi(i,j,k) = (sin(((double)i)/32 * 8*M_PI) + sin(((double)j)/32 *
8*M_PI) + sin(((double)k)/32 * 8*M_PI))/sqrt((double) 32*32*32);
                 }
     } }
      fftw_execute(plan_Phi3DForward);
               
      double mode_x, mode_y, mode_z; 
       // Get first 8 modes
        for(int m = 0; m < 8 ; m++) {
                    mode_x = 0.e0; mode_y = 0.e0; mode_z = 0.e0;
                   
         // Iterate all elements to sum up
         for(int i = 0; i <= 31; i++) {
            for(int j = 0; j <= 31; j++) {
                for(int k = 0; k <= 16; k++) {
                     // frequency symmetry around
                     mode_x += pow2(abs(phi_k(m, j, k)))/(8.f * M_PI) +
((m > 0) ? (pow2(abs(phi_k(32 -m, j, k))))/(8.f * M_PI) : 0.e0); 
                     mode_y += pow2(abs(phi_k(i, m,k)))/(8.f * M_PI) + ((m
> 0) ? (pow2(abs(phi_k(i, 32 -m, k))))/(8.f * M_PI) : 0.e0);
//complex conjugates for i > N/2 (factor 2.e0) mode_z += ((i > 0) ? 2.e0 : 1.e0) *(pow2(abs(phi_k(i, j, m))))/(8.f * M_PI); } } } std::cout << m << " " << mode_x << " " << mode_y << " " << mode_z << std::endl; } fftw_free(plan_Phi3DForward); } ========= the mode power should be equal for each direction but instead I get following result 0 31291.1 31291.1 33246.8 1 5.74174e-29 5.33214e-29 6.44005e-30 2 1.11428e-28 6.92429e-29 2.7889e-29 3 1.21869e-27 1.13776e-27 2.66253e-28 4 20860.8 20860.8 5541.14 5 7.43371e-28 5.29499e-28 1.2743e-28 6 4.48504e-28 4.74109e-28 1.38245e-28 7 8.68525e-29 1.65377e-28 3.47825e-29 for mode 0 and 4 the discrepancy can not be explained by rounding errors etc., so I guess I calculate the mode powers wrong, but I do not see why/where ? Any hints ? Thanks for reading & helping Paul
On 5 Feb, 17:03, "Mosby" <pphilscher...@gmail.com> wrote:
> Dear all, > > for quite a while I struggle with the following problem : > Having a 3-D Array of equally data in X,Y & Z , I want to calculate the > power spectra.
...
> I &#2013266080;calculate the mode powers wrong, but I do not see why/where ? Any hints > ?
I don't know what you expect to see, but be aware that the symmetry relations in ND, N > 1, are different from the simple 1D case. The symmetry is around origo. If you expect to see conjugate symmetry mirrored around the 0 axis in each transformed dimension, you will become confused. Try some simple experiments in 2D first, so you understand the effect before you move to 3D. If you have a graphics package, use it. The effect is almost counter intuitive if you only have played with 1D DFTs. Rune