Prev: Should the normal SINE and Zero padded SINE curves be the same for fft output?
Next: Should the normal SINE and Zero padded SINE curves be the samefor fft output?
From: Mosby on 5 Feb 2010 11:03
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



From: Rune Allnor on 5 Feb 2010 13:12
On 5 Feb, 17:03, "Mosby" <pphilscher...(a)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  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
 | 
Pages: 1
Prev: Should the normal SINE and Zero padded SINE curves be the same for fft output?
Next: Should the normal SINE and Zero padded SINE curves be the samefor fft output?