Path: utzoo!attcan!uunet!lll-winken!ncis!helios.ee.lbl.gov!nosc!ucsd!rutgers!mcnc!thorin!coggins!coggins
From: coggins@coggins.cs.unc.edu (Dr. James Coggins)
Newsgroups: comp.lang.c++
Subject: Re: Fast Fourier Transformer in C++
Message-ID: <6250@thorin.cs.unc.edu>
Date: 19 Jan 89 15:52:08 GMT
References: <518@cs-spool.calgary.UUCP> <6217@thorin.cs.unc.edu>
Sender: news@thorin.cs.unc.edu
Reply-To: coggins@coggins.UUCP (Dr. James Coggins)
Organization: University Of North Carolina, Chapel Hill
Lines: 416


This message contains the source code of my FFT server. It will not
compile as is without the whole rest of my library, which will be
released in March for a nominal fee (read "we're not trying to make
money on this, but we will ask you to pay for our distribution costs")
through our department's Softlab (details later).  "Without warranty"
doesn't begin to cover the disclaimer on this code.  It is posted
mainly to illustrate the technique of "process encapsulation" or
"enzymes" I mentioned in my previous article (which generated lots of
e-mail responses - You asked for it... you're gonna get it). It also
illustrates what a class looks like when constructed according to the
scheme described in my series of articles on Managing C++ Libraries
published on the Net last Nov-Dec. Greg Bollella and I are submitting
a paper on that topic to SIGPLAN Notices with some enhancements to the
techniques. 

>Some people think that a class has to be based on some data item.
>They forget that algorithm is a noun.  An algorithm that is of
>interest as an object of study itself can form a nice class. I call
>these process encapsulations, or enzymes.  My FFT server is a prime
>example, but there are many others.

To use the server, declare the server as an object with the image size
as arguments to the constructor:

    fft_server fft(256,256)

That's right, you need a different server for each image size/shape.
This server is supposed to work for images in which the row and column
sizes are powers of 2 but not necessarily the same power of 2. 

Invoke the server (with images im1, im2, im3) as follows:

    im2 = fft.forward(im1);
    im3 = fft.inverse(im2);

Usually you want to modify im2 by multiplying by a filter before
inverting.  As this is written, im3==im1 except that im1 could have
any storage format but im3 is guaranteed to be a COMPLEX image. 

One final note. When I got the big e-mail response asking to see this
code, I thought first "How nice! People are interested in seeing this
code." Then I realized that posting code would expose it (and its
undoubtedly numerous flaws) to all of you thousands of C/C++ hackers
(and you know who you are!) who could have written this more
"efficiently" or more "portably" or whatever, and might be capable of
hacking this to bits.  The last time I felt terror like that descend
on me was when I proposed (but she did say "yes").  It is time that
somebody broke through this terror barrier, and hey, I've got nothing
more interesting to do this morning. 

Seriously, I think we need a LOT more code examples posted to this
newsgroup.  Not just demonstrations of errors, but real design
examples. Reading code is not, perhaps, as viscerally involving as
writing it, but reading code can be a lot more efficient for learning
how to code than reworking design decisions that others have already
evaluated.  Newton said that he saw farther because he stood on the
shoulders of giants.  We computer people tend to stand on each other's
toes. [see also Gerald Weinberg, The Psychology of Computer
Programming].  I'll answer questions about this code as time permits.
I'll be happy to receive super-optimized versions of this server, too,
since I use it a lot. 

---------------------------------------------------------------------
Dr. James M. Coggins          coggins@cs.unc.edu
Computer Science Department   
UNC-Chapel Hill               Sticking my neck out for a good cause...
Chapel Hill, NC 27514-3175    ...again.
---------------------------------------------------------------------
 

// DEPENDENCIES FOR CLASS FFT_SERVER

#define D_FFT_SERVER
#define D_IMAGE
#define D_COMPLEX
#define D_SUBSCRIPT

#ifndef D_PRELUDE
#include "../../coolprelude.h"
#endif


// DEFINITION OF CLASS FFT_SERVER

class fft_server {

   int  axes,nrow,ncol,nbitsr,nbitsc,masklo,maskhi;
   int *bitrevr, *bitrevc;
   complex *w1,*w2;
   int *offset1,*offset2,*tr1,*tr2;

   void server_error(char*);
   void make_server(int,int);
   void fft(image);
   void fft_atno(image);
   void fft_shuffle(image);

public:

// CONSTRUCTORS FOR FFT_SERVER

   fft_server(int);
   fft_server(int,int);
   ~fft_server();
                                       
// MESSAGES FOR FFT_SERVER

   int dim()              {return axes;}
   int size()             {return nrow*ncol;}
   image forward(image);
   image inverse(image);
};

#include "fft_server.d"

fft_server::~fft_server()
{
   delete bitrevc;
   delete bitrevr;
   delete w1;
   delete w2;
   delete offset1;
   delete offset2;
}

#include "fft_server.d"

fft_server::fft_server(int s)
{
   self.make_server(s,1);
}

#include "fft_server.d"

fft_server::fft_server(int s1,int s2)
{
   self.make_server(s1,s2);
}

#include "fft_server.d"

void fft_server::fft(image im)
{
   self.fft_shuffle(im);
   self.fft_atno(im);
}

#include "fft_server.d"

void fft_server::fft_atno(image im)
{
   int npoint,nhalf,nbutter,exponent,i;
   complex temp,Wn;
   complex* buffer;
   buffer = (complex*)im.storage_addr();
   int siz=self.size();

// Perform row transforms

   npoint=1;
   int o=0;
   int p=0;
   while(npoint<ncol)
   {
      npoint*=2;
      nhalf=npoint/2;
      nbutter=siz/npoint;
      for(exponent = 0; exponent < nhalf; exponent++)
      {
         Wn = w1[p];
         p++;
         for(i = 1 ; i <= nbutter ; i++)
         {
            temp = Wn*buffer[offset2[o]];
            buffer[offset2[o]] = buffer[offset1[o]]-temp;
            buffer[offset1[o]] += temp;
            o++;
         }
      }
   }

// Perform column transforms

   npoint=1;
   p=o=0;
   while(npoint<nrow)
   {
      npoint*=2;
      nhalf=npoint/2;
      nbutter=siz/npoint;
      for(exponent = 0; exponent < nhalf; exponent++)
      {
         Wn = w2[p];
         p++;
         for(i = 1 ; i <= nbutter ; i++)
         {
            temp = Wn*buffer[tr2[o]];
            buffer[tr2[o]] = buffer[tr1[o]]-temp;
            buffer[tr1[o]] += temp;
            o++;
         }
      }
   }
}

#include "fft_server.d"

void fft_server::fft_shuffle(image im)
{

// perform 2-D bit-reverse shuffle

   int siz=self.size();
   complex* buffer = (complex*)im.storage_addr();
   complex val;
   int i,r,j,forward,reverse;

   for (i=0;i<nrow;i++)
   {
      r=i<<nbitsc;
      for (j=0;j<ncol;j++)
      {
         forward=r|j;
         reverse=(bitrevc[j]<<nbitsr)|bitrevr[i];
         if (forward<reverse)
         {
               val=buffer[forward];
               buffer[forward]=buffer[reverse];
               buffer[reverse]=val;
         }
      }
   }
}


#include "fft_server.d"

image fft_server::forward(image im)
{
   image work = im.convert(COMPLEX);
   fft(work);
   return work;
}

#include "fft_server.d"

image fft_server::inverse(image im)
{
   image work = im.convert(COMPLEX);
   subscript s(0,0,work.shape());
   for(s.init(); s.test(); s++)
   {
      work.set(s.index(),work.get_c(s.index()).conj());
   }
   fft(work);
   double t = 1.0/(double)(nrow*ncol);
   work.scale(0.0,t,0.0);
   return work;
}

#include "fft_server.d"

void fft_server::make_server(int naxis1, int naxis2)
{
   nrow=naxis1;
   ncol=naxis2;

// Determine the number of bits in the number of rows

   int num;
   num=nrow;
   nbitsr=nbitsc=0;
   while(num>1)
   {
      nbitsr++;
      num=num>>1;
   }
   num=ncol;
   while(num>1)
   {
      nbitsc++;
      num=num>>1;
   }
   masklo = nrow-1;
   maskhi = ((nrow*ncol)-1) & (~masklo);

//  This bit computes the offsets and
//  computes the Wn^e complex coefficients where
//  Wn = exp(i2pi/N) = the nth complex root of unity,
//  and e goes from 0 to N-1.

   double invar= -6.28318530717959/ncol;
   complex temp;
   axes=(naxis1==1) ? 1 : 2;
 
   int s = (int)(nrow*ncol*nbitsc/2);
   offset1 = new int[s];
   offset2 = new int[s];
   
   s = (int)(nrow*ncol*nbitsr/2);
   tr1 = new int[s];
   tr2 = new int[s];

   w1=new complex[ncol-1];
   w2=new complex[nrow-1];
 
   complex *wtemp;
   if (nrow>ncol)
      wtemp=new complex[nrow];
   else
      wtemp=new complex[ncol];
 
   for (int ex=0; ex<ncol; ex++)
   {
      temp.imag()=invar*(double)ex;
      wtemp[ex]=temp.exp();
   }

// Go through the FFT loop storing the offsets used in
// each butterfly of the row transforms
 
   int nfactor;
   int npoint=1;
   int nhalf,nbutter;
   int o=0;
   int p=0;
   while (npoint<ncol)
   {
      npoint*=2;
      nfactor=ncol/npoint;
      nhalf=npoint/2;
      nbutter=nrow*ncol/npoint;
      for (int e=0;e<nhalf;e++)
      {
         offset1[o]=e;
         offset2[o]=offset1[o]+nhalf;
         w1[p]=wtemp[e*nfactor];
         o++;
         p++;
         for (int i=1;i<nbutter;i++)
         {
            offset1[o]=offset1[o-1]+npoint;
            offset2[o]=offset2[o-1]+npoint;
            o++;
         }
      }
   }
   invar= -6.28318530717959/nrow;
   for (ex=0; ex<nrow; ex++)
   {
      temp.imag()=invar*(double)ex;
      wtemp[ex]=temp.exp();
   }

// Go through the FFT loop storing the offsets used in
// each butterfly of the row transforms

   o=0;
   p=0;
   npoint=1;
   while (npoint<nrow)
   {
      npoint*=2;
      nfactor=nrow/npoint;
      nhalf=npoint/2;
      nbutter=nrow*ncol/npoint;
      for (int e=0;e<nhalf;e++)
      {
         int temp1=e;
         int temp2=temp1+nhalf;
         tr1[o] = (((temp1&maskhi)>>nbitsr) | ((temp1&masklo)<<nbitsc));
         tr2[o] = (((temp2&maskhi)>>nbitsr) | ((temp2&masklo)<<nbitsc));
         w2[p]=wtemp[e*nfactor];
         o++;
         p++;
         for (int i=1;i<nbutter;i++)
         {
            temp1+=npoint;
            temp2+=npoint;
            tr1[o] = (((temp1&maskhi)>>nbitsr) | ((temp1&masklo)<<nbitsc));
            tr2[o] = (((temp2&maskhi)>>nbitsr) | ((temp2&masklo)<<nbitsc));
            o++;
         }
      }
   }

// This section computes the bit-reversal table

   bitrevc = new int[ncol];
   bitrevr = new int[nrow];
   int reverse,number,bits,control;

   for(number=0; number<nrow; number++)
   {
      for(reverse=0,control=nrow,bits=number; control>1;
          control=control>>1,bits=bits>>1)
         reverse = (reverse<<1) | (bits&1);
      bitrevr[number]=reverse;
   }
   for(number=0; number<ncol; number++)
   {
      for(reverse=0,control=ncol,bits=number; control>1;
          control=control>>1,bits=bits>>1)
         reverse = (reverse<<1) | (bits&1);
      bitrevc[number]=reverse;
   }
}

#include "fft_server.d"

void fft_server:: server_error(char* s)
{
   cerr << "FFT_SERVER::" << s << "\n";
   exit(2);
}