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 &#4294967295;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

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