2D Fourier Transformation in C

Question

I implemented 2D DFT and IDFT using equation from this site http://homepages.inf.ed.ac.uk/rbf/HIPR2/fourier.htm I think these are correct and nicely explained. Implementation looks like that:

    for(int i=0;i<inImage.width;i++)
    {
     for(int j=0;j<inImage.height;j++)
     {
      float ak=0; 
      float bk=0;
          for(int ii=0;ii<inImage.width;ii++)
           {
             for(int jj=0;jj<inImage.height;jj++)
              {

               float x=-2.0*PI*i*ii/(float)inImage.width;
               float y=-2.0*PI*j*jj/(float)inImage.height;
            // ak+=inImage.pixels[i][j]*(cos(x)*cos(y)-sin(x)*sin(y));
           //  bk+=inImage.pixels[i][j]*(sin(x)*cos(y)+sin(y)*cos(x));
               ak+=inImage.pixels[i][j]*cos(x+y);
               bk+=inImage.pixels[i][j]*1.0*sin(x+y);
             }
           }
     DFTImageRE.pixels[i][j]=ak;
     DFTImageIM.pixels[i][j]=bk;
       }
     }

The frequency domain (sqrt(ak * ak+bk * bk)) doesnt look as it should, and the image reconstruction (ignoring the imaginary parts) doesnt make anything near the original image. What is more pixel at [0][0] have extremely high value and no pixels range from 0 to 255 as the original one. What am i doing wrong?

Extra information:

inImage and DFTImages are just struct from which oridinary *.pgm image are construted, saving and loading images works,
i cant use any classes (like imaginary numbers) because this implementation will be on GPU side,

Thanks

@squeamishossifrage: Why should the compiler warn about conversion to a type with higher precision? Also `-Wall` does not include (all) conversion warnings. — too honest for this site, Jul 23 '16 at 14:20
Ok if you’re going to use CUDA, please for the love of all the gods, *just use [CUFFT](https://developer.nvidia.com/cufft)*! (Are you using something other than CUDA, like OpenCL?) FFT is one of the most optimized algorithms in all of computing. If your intention is to learn, then by all means implement it yourself. But if you want to get work done, please use an established library. KissFFT is a popular ANSI C lib but not the fastest. FFTW is always among the fastest open source implementations. Neither of these are that helpful in porting to GPU—use CUFFT or similar! — Ahmed Fasih, Jul 24 '16 at 00:56
@Ahmed Fasih I wil be using compute shaders to compute DFTs. Most for learning purposes so i want to implement it myself. Do you see in my implementation what could gone wrong? — Kedriik, Jul 24 '16 at 09:38
Because multiplying two complex numbers `(a + j * b) * (c + j * d)` ≠ `(a * c) + j * (b * d)` (where `j = sqrt(-1)`). Instead, `(a + j * b) * (c + j * d) = (a*c - b*d) + j * (b*c + a*d)`. Your `ak+=inImage.pixels[i][j]*cos(x+y);` is doing multiplication the first (wrong) way. If your input is all-real (image intensity), then you still have `(a + j * 0) * (c + j * d) = a*c + j * (a*d)`. — Ahmed Fasih, Jul 25 '16 at 03:05
Yes, thats why i have image for real parts (DFTImageRE) and image for imaginary parts (DFTImageIM). — Kedriik, Jul 25 '16 at 08:56

score 1 · Answer 1 · answered Aug 30 '16 at 20:15

I found a solution for my problem. It was just indexing problem. Use ii and jj in sum to find a Fourier Transform

   for(int i=0;i<inImage.width;i++)
{
 for(int j=0;j<inImage.height;j++)
 {
  float ak=0; 
  float bk=0;
      for(int ii=0;ii<inImage.width;ii++)
       {
         for(int jj=0;jj<inImage.height;jj++)
          {

           float x=-2.0*PI*i*ii/(float)inImage.width;
           float y=-2.0*PI*j*jj/(float)inImage.height;
           ak+=inImage.pixels[ii][jj]*cos(x+y);
           bk+=inImage.pixels[ii][jj]*1.0*sin(x+y);
         }
       }
 DFTImageRE.pixels[i][j]=ak;
 DFTImageIM.pixels[i][j]=bk;
   }
 }

2D Fourier Transformation in C

1 Answers1

Linked