POV-Ray : Newsgroups : povray.binaries.images : 2D Fourier Transform Server Time: 27 May 2017 08:01:51 GMT
  2D Fourier Transform (Message 1 to 3 of 3)  
From: Bald Eagle
Subject: 2D Fourier Transform
Date: 6 Mar 2017 13:10:01
Message: <web.58bd5ea376ca50b7c437ac910@news.povray.org>
After spending more than a few hours over the weekend pulling my hair out,  I've
managed to get some barely working code for a 2D transform.

Paul Bourke's source code uses "bitwise right shift" - and I don't see any
equivalent to that currently available in SDL.
Am I wrong in thinking that a bitwise right shift is just a numerical division
by 2?

I had some problems with variables in my loops - is there a limit to how many
times a variable name can be used in nested loops / includes / macros / etc?
I kept getting a "cannot assign uninitialized identifier" error when using
#local, and then "cannot pass uninitialized identifier to LValue (whatever that
is)" when using #declare.

I (sort of) fixed it - but I only needed to recode about half of the instances
of this to get the full code to work - the latter half using #local works just
fine.
Very strange.
I will post code later to illustrate.

Working on this raised many questions, that I will ask later, perhaps in
(an)other thread(s).


Post a reply to this message


Attachments:
Download '2d_fft_test_2.png' (239 KB)

Preview of image '2d_fft_test_2.png'
2d_fft_test_2.png


 

From: clipka
Subject: Re: 2D Fourier Transform
Date: 6 Mar 2017 18:17:20
Message: <58bda7b0$1@news.povray.org>
Am 06.03.2017 um 14:05 schrieb Bald Eagle:

> Paul Bourke's source code uses "bitwise right shift" - and I don't see any
> equivalent to that currently available in SDL.
> Am I wrong in thinking that a bitwise right shift is just a numerical division
> by 2?

Provided you're talking about /integer/ division by 2 (i.e. division
followed by rounding), you're perfectly right.


> I had some problems with variables in my loops - is there a limit to how many
> times a variable name can be used in nested loops / includes / macros / etc?

None that I'd be aware of.


> I kept getting a "cannot assign uninitialized identifier" error when using
> #local, and then "cannot pass uninitialized identifier to LValue (whatever that
> is)" when using #declare.

"LValue" is a moniker for "a thing you can assign values to" --
shorthand for "left-hand value", alluding to the left-hand side of an
assignment.


> I (sort of) fixed it - but I only needed to recode about half of the instances
> of this to get the full code to work - the latter half using #local works just
> fine.
> Very strange.
> I will post code later to illustrate.

Please do; I can only think that you were doing something wrong without
noticing.


Post a reply to this message

From: Bald Eagle
Subject: Re: 2D Fourier Transform
Date: 8 Mar 2017 12:55:00
Message: <web.58bffecbe7e104a8c437ac910@news.povray.org>
clipka <ano### [at] anonymousorg> wrote:

> > I (sort of) fixed it - but I only needed to recode about half of the instances
> > of this to get the full code to work - the latter half using #local works just
> > fine.
> > Very strange.
> > I will post code later to illustrate.
>
> Please do; I can only think that you were doing something wrong without
> noticing.

Here's some of the older code that I was in the middle of debugging when things
started getting weird.

###########################################################################

// Adapted to POV-Ray SDL from http://paulbourke.net/miscellaneous/dft/
#version 3.7;
global_settings {assumed_gamma 1.0}

#include "colors.inc"
//#include "functions.inc"
//#include "rand.inc"
#include "textures.inc"
#include "transforms.inc"

#declare tau = 2*pi;

default {
 texture {
  pigment {Red}
 }
}

// create a 2D data array
#declare Data2D = array [64][64][2];

#declare YPoints = dimension_size(Data2D, 1);
#declare XPoints = dimension_size(Data2D, 2);
#declare YCycles = tau/YPoints;
#declare XCycles = tau/XPoints;
#declare YAA =  2*YCycles;
#declare YBB =  4*YCycles;
#declare YCC =  6*YCycles;
#declare YDD =  8*YCycles;
#declare YEE = 10*YCycles;
#declare XAA =  1*XCycles;
#declare XBB =  3*XCycles;
#declare XCC =  5*XCycles;
#declare XDD =  7*XCycles;
#declare XEE =  9*XCycles;

#for (Y, 0, dimension_size(Data2D, 1)-1 )
 #for (X, 0, dimension_size(Data2D, 2)-1 )
  #declare Data2D[Y][X][0] =  sin(YAA*Y) + sin(XAA*X) +
        sin(YBB*Y) + sin(XBB*X) +
        sin(YCC*Y) + sin(XCC*X) +
        sin(YDD*Y) + sin(XDD*X) +
        sin(YEE*Y) + sin(XEE*X);
  #declare Data2D[Y][X][1] = 0; // imaginary portion
 #end // end for X
#end // end for Y

//##########################################################################

#macro PowerOfTwo (N) //, optional M, optional TwoPm)
 //TwoPm ("Two to the Power of M") = 2^M, which is <= N
 //#ifndef(local.M) #local M = 0; #end
 //#ifndef(local.TwoPm) #local TwoPm = 1; #end
 #if (N <= 1)
  #declare M = 0;
  #declare TwoPm = 1;
  false
 #end // end if

 #declare M = 1;
 #declare TwoPm = 2;
 #while (2*TwoPm <= N)
  #declare M = M +1;
  #declare TwoPm = TwoPm * 2;
 #end // end while

 #if (TwoPm != N)
  false;
 #else
  true;
 #end // end if
#end // end macro PowerofTwo

//##########################################################################

#macro FFT (Direction, MM, _X, _Y, Verbose)
 // MM is calculated in the "PowerofTwo" macro
 #local NN = 1;
 #for (i, 0, MM-1)
  #local NN = NN*2;
  #if (Verbose)
   #debug concat ("NN =", str (NN, 3, 0), "\n")
  #end // end if Verbose
 #end
 // bit reversal
 //i2 = NN >>1;  // bitwise right shift
 #local i2 = NN/2;
 #declare j = 0;
 #for (i, 0, NN-1)
  #if (i < j)
   #if (Verbose)
    #debug concat ("i =", str (i, 3, 0), "     j =", str (j, 3, 0), "\n")
   #end // end if Verbose
   #local TX = _X[i];
   #local TY = _Y[i];
   #local _X[i] = _X[j];
   #local _Y[i] = _Y[j];
   #local _X[j] = TX;
   #local _Y[j] = TY;
  #end // end if
  #local K = i2;
  #while (K < j)
   #local j = j-K;
   // K >>= 1;  // assignment by bitwise right shift
   #local K = K/2;
  #end // end while
  #local j = j+K;
 #end // end for i

 // compute FFT
 #local C1 = -1;
 #local C2 = 0;
 #local L2 = 1;
 #for (L, 0, MM)
  #local L1 = L2;
  // L2 <<= 1;
  #local U1 = 1;
  #local U2 = 0;
  #for (j, 0, L1)
   #for (i, j, NN, L2)
    #local I1 = i + L1;
    #local T1 = U1 * _X[I1] - U2 * _Y[I1];
    #local T1 = U1 * _Y[I1] + U2 * _X[I1];
    #declare _X[I1] = _X[i] - T1;
    #declare _Y[I1] = _Y[i] - T2;
    #declare _X[i] = _X[i] + T1;
    #declare _Y[i] = _Y[i] + T2;
   #end // end for i
   #local  Z = U1 * C1 - U2 * C2;
   #local U2 = U1 * C2 + U2 * C1;
   #local U1 = Z;
  #end // end for j
  #local C2 = sqrt ((1-C1)/2);
  #if (Direction = 1)
   #local C2 = -C2;
  #end // end if
  #local C1 = sqrt ((1+C1)/2);
 #end // end for L

 // Scaling for forward transform
 #if (Direction = 1)
  #for (i, 0, NN)
   #declare _X[i] = _X[i] / NN;
   #declare _Y[i] = _Y[i] / NN;
  #end // end if
 #end // end if

 true

#end // end macro FFT

#macro FFT2D (Data2DArray, Direction, Verbose)
 //Perform an in-place 2D FFT given a complex 2D array
 //The value for Direction is 1 for forward, -1 for reverse
 //The size of the array is NX x NY
 //Return false if the dimensions are not powers of 2
 #local NY = dimension_size(Data2DArray, 1);
 #local NX = dimension_size(Data2DArray, 2);
 #if (Verbose)
  #debug "\n\n######################################################\n"
  #debug concat ("Starting 2D Fourier Transform of", str(NX, 3, 0), " x",
str(NY, 3, 0), " sample matrix \n")
  #debug concat ("Transforming", str(NX, 3, 0), " Rows... \n")
 #end // end if Verbose
 #local P2 = PowerOfTwo(NX) //, M, TwoPm);
 #if (!P2 | TwoPm != NX)
  false;
 #end
 #local Row_Real = array [NX];
 #local Row_Imaginary = array [NX];
 #local C_Real = array [NY][NX];
  #for (j, 0, NX-1)
   #for (i, 0, NY-1)
    #local C_Real[i][j] = Data2DArray[i][j][0];
   #end // end for j
  #end // end for i
 #local C_Imaginary = array [NY][NX];
  #for (j, 0, NX-1)
   #for (i, 0, NY-1)
    #local C_Imaginary[i][j] = Data2DArray[i][j][1];
   #end // end for j
  #end // end for i
 #for (j, 0, NY-1)
  #for (i, 0, NX-1)
   #local Row_Real[i] = C_Real[i][j];
   #local Row_Imaginary[i] = C_Imaginary[i][j];
  #end // end for i
  FFT(Direction, M, Row_Real, Row_Imaginary, yes)
  #for (i, 0, NX-1)
   #local C_Real[i][j] = Row_Real[i];
   #local C_Imaginary[i][j] = Row_Imaginary[i];
  #end // end for i
 #end // end for j

 #if (Verbose)
  #debug concat ("Transforming ", str(NY, 3, 0), " Columns...\n")
 #end // end if Verbose
 #local P2 = PowerOfTwo(NY)
 #if (!P2 | TwoPm != NY)
  false
 #end
 #declare Col_Real = array [NY];
 #declare Col_Imaginary = array [NY];
 //#declare C_Real = array [NY][NX];
 //#declare C_Imaginary = array [NY][NX];
 #for (i, 0, NX-1)
  #for (j, 0, NY-1)
   #declare Col_Real[j] = C_Real[i][j];
   #declare Col_Imaginary[j] = C_Imaginary[i][j];
  #end // end for i
  FFT(Direction, M, Col_Real, Col_Imaginary, yes)
  #for (j, 0, NY-1)
   #declare C_Real[i][j] = Col_Real[j];
   #declare C_Imaginary[i][j] = Col_Imaginary[j];
  #end // end for j
 #end // end for i
 #if (Verbose)
  #debug "2D Fourier Transform complete."
 #end // end if Verbose
 true
#end // end macro FFT2D

#declare FFTArray = FFT2D (Data2D, 1, yes)




/*
int FFT2D(COMPLEX **c, int nx, int ny, int dir)
{
   int i,j;
   int m,twopm;
   double *real,*imag;

   // Transform the rows
   real = (double *)malloc(nx * sizeof(double));
   imag = (double *)malloc(nx * sizeof(double));
   if (real == NULL || imag == NULL)
      return(FALSE);
   if (!Powerof2(nx,&m,&twopm) || twopm != nx)
      return(FALSE);
   for (j=0;j<ny;j++) {
      for (i=0;i<nx;i++) {
         real[i] = c[i][j].real;
         imag[i] = c[i][j].imag;
      }
      FFT(dir,m,real,imag);
      for (i=0;i<nx;i++) {
         c[i][j].real = real[i];
         c[i][j].imag = imag[i];
      }
   }
   free(real);
   free(imag);

   // Transform the columns
   real = (double *)malloc(ny * sizeof(double));
   imag = (double *)malloc(ny * sizeof(double));
   if (real == NULL || imag == NULL)
      return(FALSE);
   if (!Powerof2(ny,&m,&twopm) || twopm != ny)
      return(FALSE);
   for (i=0;i<nx;i++) {
      for (j=0;j<ny;j++) {
         real[j] = c[i][j].real;
         imag[j] = c[i][j].imag;
      }
      FFT(dir,m,real,imag);
      for (j=0;j<ny;j++) {
         c[i][j].real = real[j];
         c[i][j].imag = imag[j];
      }
   }
   free(real);
   free(imag);

   return(TRUE);
}
*/

/*-------------------------------------------------------------------------
   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

     Formula: forward
                  N-1
                  ---
              1   \          - j k 2 pi n / N
      X(n) = ---   >   x(k) e                    = forward transform
              N   /                                n=0..N-1
                  ---
                  k=0

      Formula: reverse
                  N-1
                  ---
                  \          j k 2 pi n / N
      X(n) =       >   x(k) e                    = forward transform
                  /                                n=0..N-1
                  ---
                  k=0
*/


/*
int FFT(int dir,int m,double *x,double *y)
{
   long nn,i,i1,j,k,i2,l,l1,l2;
   double c1,c2,tx,ty,t1,t2,u1,u2,z;

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

   // Do the bit reversal
   i2 = nn >> 1;
   j = 0;
   for (i=0;i<nn-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;
      } // end if
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      } // end while
      j += k;
   } // end for i

   // 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<nn;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;
         } // end for i
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      } // end for j
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1)
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   } // end for L

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

   return(TRUE);
}
*/

/*-------------------------------------------------------------------------
   Calculate the closest but lower power of two of a number
   twopm = 2**m <= n
   Return TRUE if 2**m == n

int Powerof2(int n,int *m,int *twopm)
{
   if (n <= 1) {
      *m = 0;
      *twopm = 1;
      return(FALSE);
   }

   *m = 1;
   *twopm = 2;
   do {
      (*m)++;
      (*twopm) *= 2;
   } while (2*(*twopm) <= n);

   if (*twopm != n)
      return(FALSE);
   else
      return(TRUE);
}
*/


Post a reply to this message

Copyright 2003-2008 Persistence of Vision Raytracer Pty. Ltd.