Thread: Fast Fourier Transforms, k values?

Hi there,

I getting very confused. I have written a simple program that compares a Gaussian in real space with a Gaussian that is reverse Fourier transformed from Fourier space.
My goal is to work out what the k values should be for each element of my Gaussian array in Fourier space.
In 1D I have the Gaussian 1/(sqrt(2*pi)*RG*RG) * (exp(-(i*i)/(2*RG*RG)) + exp(-(i-(N-1))*(i-(N-1))/(2*RG*RG)),

where RG is my Gaussian filter size and i is the x value which runs from 0N-1).

In k-space, I have a Gaussian G(k) = exp(-kx*kx*RG*RG/2), where kx runs from (2*pi/N) * [0N/2 -1) -N/2 -(N/2-1):-1].
when I reverse FFT and plot the output they look the same.

However, when I move to 2D this goes wrong.

So, as a check of the kx values, I FFT'd the real Gaussian and tried to extract the kx from there, i.e kx = sqrt(2*log(FFT(G(r)*N))/(RG*RG), and did not get the same kx values as I used in the 'raw' Gaussian G(kx).

Can anyone help me with this?

Thank you!

Err the frown icon is meant to be a colon..

Post your code... you'll get a lot more help than posting a jiant string of math at us.

Code:

#include "header.h"
#include<stdio.h>
#include<fftw3.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include <malloc.h>
#include <fstream>
#include <iomanip>
#include<complex.h>


const int NG = 10;
const double NG_i = NG*((NG/2) +1);
const double NGMAX = NG*NG;
const double imag = sqrt(-1);
double SQ(double);
double CU(double);

int main() {

/*~~~~~~~~~~~~~~~First use the real Gaussian and FFT to convolve~~~~~~~~~~*/
 
  double *gauss = (double*)malloc(sizeof(double)*NGMAX);
  for (int i =0;i<NG;i++) {
    for (int j=0;j<NG;j++) {
      int m = j + i*NG;
      int rs = i*i + j*j;
      int rsd = (i)*(i) + (j)*(j);
      int rsg2 = (i*i) + (j-(NG-1))*(j-(NG-1));
      int rsg3 = (i-(NG-1))*(i-(NG-1)) + (j-(NG-1))*(j-(NG-1));
      int rsg4 = (i-(NG-1))*(i-(NG-1)) + j*j;
      gauss[m] = (double)(1/((2*pi)))*(exp(-0.5*rs) 
                 + exp(-0.5*rsg2)
                 + exp(-0.5*rsg3)
                 + exp(-0.5*rsg4));
    }
  }
  /*save this data*/
  char Fout_name[200];
  sprintf(Fout_name,"gauss1.dat");
  FILE *Fout;
  if((Fout=fopen(Fout_name,"w"))==NULL) printf("conv 1 data file not opened");
  for(int idd=0; idd<NGMAX; idd++) {
      fprintf(Fout,"%lf\n",gauss[idd]);

  }
  fclose(Fout);
  for (int i=0;i<5;i++) {
    printf("gauss1 [%d] = %lf\n",i,gauss[i]);
  }

  /*fft the gaussian*/
  fftw_complex *fft_result;
  fftw_plan p; 
  fft_result  = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * NGMAX);
  p  = fftw_plan_dft_r2c_2d(NG,NG, gauss, fft_result, FFTW_ESTIMATE );
  fftw_execute( p );
 
  fftw_destroy_plan(p);
  
  printf("1\n");

  /*then reverse the transform to get the real gaussian back*/
  double  *gauss1b;
  fftw_plan p1; 
  gauss1b = ( double* )malloc( sizeof( double ) * NGMAX);
  p1  = fftw_plan_dft_c2r_2d(NG,NG,fft_result, gauss1b, FFTW_ESTIMATE );
  fftw_execute( p1);

  printf("3\n");

  /*free memory*/
  fftw_free(fft_result);
  fftw_destroy_plan(p1);

  /*save this data*/
  char Fout_name2[200];
  sprintf(Fout_name2,"gauss1b.dat");
  FILE *Fout2;
  if((Fout2=fopen(Fout_name2,"w"))==NULL) printf("conv 1 data file not opened");
  for(int idd=0; idd<NGMAX; idd++) {
      fprintf(Fout2,"%lf\n",gauss1b[idd]/NGMAX);
  }
  fclose(Fout2);
  for (int i=0;i<5;i++) {
    printf("gauss1 trans [%d] = %lf\n",i,gauss1b[i]/NGMAX);
  }
  printf("4\n");

/*~~~~~~~~~~~~~~~Now use the Gaussian(k) as a comparison~~~~~~~~~~~~~~~~~~~~~*/

  /*compute the 'raw' Fourier using k values as a check of k*/
  double kx,ky,fy;
  fftw_complex *gauss2; 
  gauss2 = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*NGMAX);
  for (int i=0;i<NG;i++) {
    for (int j=0;j<NG;j++) {
      if (i<NG/2) kx =(double)i*2*pi/((double)NG);
      else kx = ((double)i-(double)NG)*2*pi/((double)NG);
      if (j<NG/2) ky =(double)j*2*pi/((double)NG);
      else ky = ((double)j-(double)NG)*2*pi/((double)NG);
      int m = j + i*(NG);
      double ks = (kx*kx) + (ky*ky);
      gauss2[m][0] = (exp(-ks/2));
      gauss2[m][1] =0.0;
      
    }
  }

  printf("5\n");


  /*then reverse the transform to get the real gaussian*/
  double  *gauss2re;
  fftw_plan p2; 
  gauss2re = ( double* )malloc( sizeof( double ) * NGMAX);
  p2 = fftw_plan_dft_c2r_2d(NG,NG,gauss2, gauss2re, FFTW_ESTIMATE );
  fftw_execute( p2);

  /*free memory*/
  fftw_free(gauss2);
  fftw_destroy_plan(p2);

  /*save this data*/
  char Fout_name3[200];
  sprintf(Fout_name3,"gauss2b.dat");
  FILE *Fout3;
  if((Fout3=fopen(Fout_name3,"w"))==NULL) printf("conv 1 data file not opened");
  for(int idd=0; idd<NGMAX; idd++) {
      fprintf(Fout3,"%lf\n",gauss2re[idd]/(double)NGMAX);
  }
  fclose(Fout3);
  for (int i=0;i<5;i++) {
    printf("Fourier Gauss [%d] = %lf\n",i,gauss2re[i]/(double)NGMAX);
  }

}
double SQ(double N) {
  double Ns = N*N;
  return Ns;
}

double CU(double N) {
  double Nc = N*N*N;
  return Nc;
}

All I am trying to do is find which values of k correspond to each element of the Fourier transformed function and I want a way of checking them so that I am sure.
I have written the Gaussian in k space as this requires k values. I wanted to check if when I reverse transformed it, it would be the same as the Gaussian in real space. When I print out the elements of each array they do not match. Thus, I am guessing that my k values are wrong.