Multiple Reverse FFTs don't work?

Hello,

I am trying to impliment a simple FFT algorithm I found on the web. 
Algorithm is <a
href="http://astronomy.swin.edu.au/~pbourke/other/dft/">here.</a><P>  

When I set up a single 1 dimensional 'row' of (any) data, I can perform
the FFT on it, and then perform the reverse FFT on it and the final output
is identical to the input.  

What I want to do is perform a 2d fft on a grayscale image.  But the first
thing I am doing (as a test) is failing.  I am simply using the algorithm
above to go through a grayscale image, line by line and FFT the line, then
store the transform in another 2D array.  

Then, I am going back through that stored array, and reverse-FFT ing each
line.  

It seems like this should return the correct initial image, since its just
repeating the 'single row' behavior.  However, the image outputted is all
wrong.  Its correct for a couple lines, then becomes noise, then nothing
at all.

Does anyone have any idea what I'm doing wrong.  This is driving me
crazy!!!

typewriter wrote:

> Hello,
> 
> I am trying to impliment a simple FFT algorithm I found on the web. 
> Algorithm is <a
> href="http://astronomy.swin.edu.au/~pbourke/other/dft/">here.</a><P>  
> 
> When I set up a single 1 dimensional 'row' of (any) data, I can perform
> the FFT on it, and then perform the reverse FFT on it and the final output
> is identical to the input.  
> 
> What I want to do is perform a 2d fft on a grayscale image.  But the first
> thing I am doing (as a test) is failing.  I am simply using the algorithm
> above to go through a grayscale image, line by line and FFT the line, then
> store the transform in another 2D array.  
> 
> Then, I am going back through that stored array, and reverse-FFT ing each
> line.  
> 
> It seems like this should return the correct initial image, since its just
> repeating the 'single row' behavior.  However, the image outputted is all
> wrong.  Its correct for a couple lines, then becomes noise, then nothing
> at all.
> 
> Does anyone have any idea what I'm doing wrong.  This is driving me
> crazy!!!
> 
> 
Your theory is correct, so there must be something wrong with your 
algorithm or your code.  I can't say what it is, though.

Two things you can try are to either find someone who owes you a favor 
and get them to go through the code, or to simplify your problem -- what 
happens if you do the FFT on the same data repeatedly?  Does it 
eventually not return the right answer?

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Posting from Google?  See http://cfaj.freeshell.org/google/

"Applied Control Theory for Embedded Systems" came out in April.
See details at http://www.wescottdesign.com/actfes/actfes.html

typewriter wrote:

> I am trying to impliment a simple FFT algorithm I found on the
> web. Algorithm is <a
> href="http://astronomy.swin.edu.au/~pbourke/other/dft/">here.</a>
> <P>

Usenet ain't a web board, whatever Google would have you think.

So you are using the code in Bourke's Appendix B?

> When I set up a single 1 dimensional 'row' of (any) data, I can
> perform the FFT on it, and then perform the reverse FFT on it
> and the final output is identical to the input.  
> 
> What I want to do is perform a 2d fft on a grayscale image.  But
> the first thing I am doing (as a test) is failing.  I am simply
> using the algorithm above to go through a grayscale image, line
> by line and FFT the line, then store the transform in another 2D
> array.  
> 
> Then, I am going back through that stored array, and reverse-FFT
> ing each line.  
> 
> It seems like this should return the correct initial image,
> since its just repeating the 'single row' behavior.  However,
> the image outputted is all wrong.  Its correct for a couple
> lines, then becomes noise, then nothing at all.

It's hard to diagnose such problems without seeing the salient parts 
of your code. However, the progressive worsening suggests to me that 
you might not be zero-padding the rows to a power-of-two length so 
that each FFT call actually works on the current row *plus* a prefix 
of the next row; or you are not encoding the input as pairs of double  
with contents (pixel[n], 0.0); or even both.


Martin

-- 
Quidquid latine scriptum sit, altum viditur.

>It's hard to diagnose such problems without seeing the salient parts 
>of your code. 
.
>you might not be zero-padding the rows to a power-of-two length 
.
>or you are not encoding the input as pairs of double  

>Martin

Thanks for your reply.

Woops sorry about the HTML.  Anyway, the data is 2^9 512 pixels/row, so
there is no zero padding nessesary.  And yes I am zeroing out the
imaginary part of the input.  The pixel data is float grayscale [0.0,1.0].
  Does FFT algorithm care about the scale of the signal?  ie [0,1.] vs.
[0,256]?  

Here is the code, it is written in the Processing interpreted language,
which I beleive follows the basic syntax of Java.

//*************************************

float data2[]=new float[512];
float im[]=new float[512*512];
float im2[]=new float[512*512];
float temp[]=new float [512];
  
int i,j;
float ii,jj;
  
void setup(){
 size(512,512);
 colorMode(RGB,1.0);
 background(0,0,0);
  for(j=0;j<512;j++){         /// create an image in im
    for(i=0;i<512;i++){   
      ii=i;jj=j;
      im[(512*j)+i]=(abs(cos(dist(ii,jj,0,0)/5.)));  // Circles      
    }
    data2[j]=0;       //zero out data2
  }
    
  
}



void draw(){         //MAIN LOOP
  noLoop();
 
  fft2d(im,im,1);    // Subroutine: in place row by row forward FFT of im
  fft2d(im,im,-1);   // Subroutine: in place row by row reverse FFT of im
  
  for(j=0;j<511;j++){  // DRAW CONTENTS OF im[]
    for(i=0;i<511;i++){
      stroke(im[j*width+i]);
      point(i,j);
    }
  }
}


void fft2d(float in[], float out[], int direction){
  int i,j;

  
  for(j=0;j<512;j++){         // For each row
    for(i=0;i<512;i++){       // Go thru the row
      temp[i]=in[j*512+i];    // Fill temp[] with row contents
    }

    FFT(direction,9,temp,data2);  // In place FFT of temp[]
    // ^---1D version, which tests seem to show as working 

    for(i=0;i<512;i++){       // Go thru the row
      out[j*512+i]=temp[i];   // Put row into correct place in out[]
    }
  }

//      for(i=0;i<512;i++){data2[i]=0;}   //Re-zero data2[]????

  
}

/*
   This computes an in-place complex-to-complex FFT 
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform 
*/
int FFT(int dir,int m,float []x,float []y)
{
   int n,i,i1,j,k,i2,l,l1,l2;
   float c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   n = 1;
   for (i=0;i<m;i++) 
      n *= 2;

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for (i=0;i<n-1;i++) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* Compute the FFT */
   c1 = -1.0; 
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0; 
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1; 
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1) 
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<n;i++) {
         x[i] /= n;
         y[i] /= n;
      }
   }
   
   return(1);
}

typewriter wrote:

> Does FFT algorithm care about the scale of the signal?  ie
> [0,1.] vs. [0,256]?  

I don't think so. The DFT is linear and you have no fixed-point 
worries.

> Here is the code, it is written in the Processing interpreted
> language, which I beleive follows the basic syntax of Java.
> 
> //*************************************
> 
> float data2[]=new float[512];
> float im[]=new float[512*512];
> float im2[]=new float[512*512];
> float temp[]=new float [512];
>   
> int i,j;
> float ii,jj;
>   
> void setup(){
>  size(512,512);
>  colorMode(RGB,1.0);
>  background(0,0,0);
>   for(j=0;j<512;j++){         /// create an image in im
>     for(i=0;i<512;i++){   
>       ii=i;jj=j;
>       im[(512*j)+i]=(abs(cos(dist(ii,jj,0,0)/5.)));  // Circles 
>           
>     }
>     data2[j]=0;       //zero out data2
>   }
>     
>   
> }
> 
> 
> 
> void draw(){         //MAIN LOOP
>   noLoop();
>  
>   fft2d(im,im,1);    // Subroutine: in place row by row forward
>   FFT of im fft2d(im,im,-1);   // Subroutine: in place row by
>   row reverse FFT of im 
>   
>   for(j=0;j<511;j++){  // DRAW CONTENTS OF im[]
>     for(i=0;i<511;i++){
>       stroke(im[j*width+i]);
>       point(i,j);
>     }
>   }
> }
> 
> 
> void fft2d(float in[], float out[], int direction){

You might consider another name as a "2D" FFT would usually be 
understood as transforming along both spatial directions.

>   int i,j;
>   
>   for(j=0;j<512;j++){         // For each row
>     for(i=0;i<512;i++){       // Go thru the row
>       temp[i]=in[j*512+i];    // Fill temp[] with row contents
>     }
> 
>     FFT(direction,9,temp,data2);  // In place FFT of temp[]
>     // ^---1D version, which tests seem to show as working 

At this point, the DFT of the j-th row is the vector temp + I*data2. 
You go on to discard the imaginary part for output, so the function 
does not actually return the full row FFT leaving you unable to 
recover the image.

Maybe what you intended was discarding half of the bins to capitalize 
on real-input Hermitian symmetry? But then the inverse transform 
would have to take that into account.

>     for(i=0;i<512;i++){       // Go thru the row
>       out[j*512+i]=temp[i];   // Put row into correct place in
>                               // out[] 
>     }
>   }
> 
> //      for(i=0;i<512;i++){data2[i]=0;}   //Re-zero data2[]????

Yes, you need to zero the imaginary part for each row. Note that this 
line is at the wrong nesting level as written.

> }


Martin

-- 
Quidquid latine scriptum sit, altum viditur.

Nice!  Thanks very much.  Something about what you wrote triggered
understanding.  I was unclear that a real series DFT created an imaginary
component that was NEEDED to later recover the space domain information. 
My fft2d func now fits its name.

One last question.  Now I can forward DFT and Reverse DFT and get the
original image back just fine.  But how do I view the freq. domain image. 
As of now, it is almost totally black.  Is this where I need to compute the
vector using both the imag. and real components?


Thanks again man, 
Kristopher 


Here's the working code:

//**************************************************88
float real[]=new float[512*512];  // image real component (pixel data)
float imag[]=new float[512*512];  // image imag component 

float rtemp[]=new float [512];    // Single row real components
float itemp[]=new float [512];    // single row imag components

  
int i,j;
float ii,jj;
  
void setup(){
 size(512,512);
 colorMode(RGB,1);
 background(0,0,0);
  for(j=0;j<512;j++){         
    for(i=0;i<512;i++){   
      ii=i;jj=j;
      real[(512*j)+i]=(abs(cos(dist(ii,jj,0,0)/5.)));        // define
pixels in real part
      imag[(512*j)+i]=0;                                    // zero out
imag part
    }
  }
}



void draw(){
  noLoop();
 
  fft2d(real,imag,1);
  fft2d(real,imag,-1);
 
  for(j=0;j<511;j++){  
    for(i=0;i<511;i++){
      stroke(real[j*width+i]);
      point(i,j);
    }
  }
}


void fft2d(float real[], float imag[], int direction){
    int i,j;
// FIRST DO ALL THE ROWS
    for(j=0;j<512;j++){
      for(i=0;i<512;i++){
        rtemp[i]=real[j*512+i];      // Load rtemp[] with row data 
        itemp[i]=imag[j*512+i];
      }

      FFT(direction,9,rtemp,itemp);  // In place DFT of temp row data 

      for(i=0;i<512;i++){
        real[j*512+i]=rtemp[i];      // store real and imaginary data in
large buffers
        imag[j*512+i]=itemp[i];
      }
    }
    for(i=0;i<512;i++){      itemp[i]=0;    }    
// NOW DO VERTICAL COLUMNS
    for(i=0;i<512;i++){
      for(j=0;j<512;j++){
        rtemp[j]=real[j*512+i];      // Load rtemp[] with col data
        itemp[j]=imag[j*512+i];
      }

      FFT(direction,9,rtemp,itemp);  // In place DFT of temp row data 

      for(j=0;j<512;j++){
        real[j*512+i]=rtemp[j];      // store real and imaginary data in
large buffers
        imag[j*512+i]=itemp[j];
      }
    }
    for(i=0;i<512;i++){      itemp[i]=0;    }    


}

/*
   This computes an in-place complex-to-complex FFT 
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform 
*/
int FFT(int dir,int m,float []x,float []y)
{
   int n,i,i1,j,k,i2,l,l1,l2;
   float c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   n = 1;
   for (i=0;i<m;i++) 
      n *= 2;

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for (i=0;i<n-1;i++) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* Compute the FFT */
   c1 = -1.0; 
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0; 
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1; 
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1) 
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<n;i++) {
         x[i] /= n;
         y[i] /= n;
      }
   }
   
   return(1);
}

typewriter wrote:

> Nice!  Thanks very much.  Something about what you wrote
> triggered understanding.  I was unclear that a real series DFT
> created an imaginary component that was NEEDED to later recover
> the space domain information. My fft2d func now fits its name.

Great!

> One last question.  Now I can forward DFT and Reverse DFT and
> get the original image back just fine.  But how do I view the
> freq. domain image. As of now, it is almost totally black.  Is
> this where I need to compute the vector using both the imag. and
> real components? 

It depends on what you're interested in whether it makes the most 
sense to display a complex buffer as real and imaginary parts or as 
magnitude and phase. In each case you can just do two plots or use 
something like brightness and color to put all the information in one 
plot. Either way, I'd suggest you scale each component by its maximum 
value to use the full display resolution, indicating the scale on 
your UI.


Martin

-- 
Quidquid latine scriptum sit, altum viditur.

typewriter wrote:

> void fft2d(float real[], float imag[], int direction){

There's no need to zero itemp any more since you load it from imag on 
each iteration anyway.

You have a lot of copying going on here. It might be worthwile to 
come up with a View abstraction that wraps an array with a 
configurable indexing scheme, allowing subarrays and non-unit 
strides, then use that facility to present chunks of real and imag to 
the FFT and compare runtimes.


Martin

-- 
Quidquid latine scriptum sit, altum viditur.

>There's no need to zero itemp any more since you load it from imag on 
>each iteration anyway.

ah yes..

>You have a lot of copying going on here. It might be worthwile to ...

Yes, very good suggestions.  I'll probably work on some of that when I 
really know whats going on with this stuff.

So far I've had a couple of observations.  First, the Reverse DFT 
doesn't seem to care whether the frequency spectrum quadrants have been 
swapped.  The image is recreated identically with or without quadrant 
swapping.  

Second, to view the power spectrum, I calculated the lenth of the vector
of
each pixel, using the real and imaginary parts.  Then applied logarithmic

scaling to push the values up into a more viewable range.  It looks nice. 

Of course, I dont touch the actual data - its still needed for the reverse

DFT.

Performing filtering is like _magic_.  Its amazing to play with.  
I'm really stoked to be at this point.  Thanks 
Martin for your help.   Soon I am hoping to move into some color 
experiments by DFTing R,G and B separately....

typewriter wrote:

> So far I've had a couple of observations.  First, the Reverse
> DFT doesn't seem to care whether the frequency spectrum
> quadrants have been swapped.

Not sure what you mean by quadrant swapping?

> Second, to view the power spectrum, I calculated the lenth of
> the vector of each pixel, using the real and imaginary parts.

Better reserve the term "pixel" for the spatial domain and call the 
DFT points something else. It think people still use "bin" in the 2D 
case.


Martin

-- 
Quidquid latine scriptum sit, altum viditur.