[comp.lang.c++] Fast Fourier Transformer in C++

coggins@coggins.cs.unc.edu (Dr. James Coggins) (01/19/89)

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);
}

samperi@marob.MASA.COM (Dominick Samperi) (01/22/89)

In article <6250@thorin.cs.unc.edu> coggins@coggins.UUCP (Dr. James Coggins) writes:
>Seriously, I think we need a LOT more code examples posted to this
>newsgroup.  Not just demonstrations of errors, but real design
>examples. 

Well, this isn't a real design example, but it addresses an important C++
issue, namely, parametrized classes, implemented via the macro package
in <vector.h>. First, it illustrates that there are some rather messy
syntactical restrictions (no blanks in the macro formal parameter list,
macro expansion inside of comments, etc.), and second, it illustrates
that there are some other problems, leading to bugs when the macro
package is used (which is the case when BUG is defined). The bug is
described in counter::counter(counter& cnt).

Any insights would be appreciated.

Dominick Samperi
samperi@acf8.nyu.edu
uunet!hombre!samperi

/* #define BUG /* Uncomment this to use <vector.h> and cause a BUG! */

#include <stream.h>

#ifdef BUG 

#include <vector.h>

declare(vector,int) ;	// NOTE: no spaces allowed between () here!
implement(vector,int) ;	// NOTE: this macro should only be expanded in ONE file.
			//       ^^^^ can't say the implXXX word here, since
			//            the macro processor doesn't understand //!

#else

class vector {
	int *v ;
	int sz ;
public:
	vector(int) ;
	vector(vector&) ;
	~vector() ;
	int size()	{ return sz ; }
	int& operator[](int) ;
	vector operator=(vector) ;
} ;

vector::vector(int siz)
{
	sz = siz ;
	v = new int[sz] ;
}

vector::vector(vector& vv)
{
	sz = vv.size() ;
	v = new int[sz] ;
	for(int i = 0 ; i < sz ; i++) v[i] = vv[i] ;
}

vector::~vector()
{
	delete v ;
}

int& vector::operator[](int i)
{
	return v[i] ;
}

vector vector::operator=(vector vv)
{
	for(int i = 0 ; i < sz ; i++) v[i] = vv[i] ;
	return *this ;
}

#endif

class counter {
#ifdef BUG
	vector(int) *low, *high, *incr, *count ;
#else
	vector *low, *high, *incr, *count ;
#endif
	int sz ;
public:
#ifdef BUG
	counter(vector(int)&, vector(int)&, vector(int)&) ;
#else
	counter(vector&, vector&, vector&) ;
#endif
	counter(counter&) ;
	~counter() ;
	int size() { return low->size() ; }
	void reset()	{ *count = *low ; }
	void print() ;
	int operator[](int i)	{ return (*count)[i] ; }
} ;

void counter::print()
{
	cerr << "Counter(low,high,incr,count): \n" ;
	for(int i = 0 ; i < sz ; i++)
	{
		cerr << i << ": " << (*low)[i] << " " << (*high)[i] << " "
		     << (*incr)[i] << " " << (*count)[i] << "\n" ;
	}
}

#ifdef BUG
counter::counter(vector(int)& lo, vector(int)& hi, vector(int)& inc)
#else
counter::counter(vector& lo, vector& hi, vector& inc)
#endif
{
	sz = lo.size() ;
#ifdef BUG
	low = new vector(int)(sz) ;
	high = new vector(int)(sz) ;
	incr = new vector(int)(sz) ;
	count = new vector(int)(sz) ;
#else
	low = new vector(sz) ;
	high = new vector(sz) ;
	incr = new vector(sz) ;
	count = new vector(sz) ;
#endif
	*low = lo ;
	*high = hi ;
	*incr = inc ;
	*count = lo ;
}

counter::counter(counter& cnt)
{

cerr << "cnt param before instantiation:\n" ;
cnt.print() ;

	sz = cnt.size() ;
#ifdef BUG

/*
*  When BUG is defined, these memory allocation statements cause cnt to be
*  corrupted, as you will see from the output generated by cnt.print() before
*  and after this code.
*/

	low = new vector(int)(sz) ;
	high = new vector(int)(sz) ;
	incr = new vector(int)(sz) ;
	count = new vector(int)(sz) ;
#else
	low = new vector(sz) ;
	high = new vector(sz) ;
	incr = new vector(sz) ;
	count = new vector(sz) ;
#endif

cerr << "cnt after instantiation:\n" ;
cnt.print() ;

	*low = *cnt.low ;
	*high = *cnt.high ;
	*incr = *cnt.incr ;
	*count = *low ;
}

counter::~counter()
{
	delete count ;
	delete incr ;
	delete high ;
	delete low ;
}

ostream& operator<<(ostream& s, counter& c)
{
	s << "[ " ;
	for(int i = 0 ; i < c.size() ; i++)
		s << form("%2d", c[i]) << " " ;
	return s << "]\n" ;
}

main()
{
#ifdef BUG
	vector(int) low(3), high(3), incr(3) ;
#else
	vector low(3), high(3), incr(3) ;
#endif

	for(int i = 0 ; i < 3 ; i++)
	{
		low[i] =  1 ;
		high[i] = 5 ;
		incr[i] = 1 ;
	}

	counter select(low, high, incr) ;

	select.print() ;

	counter x = select ;
}

jima@hplsla.HP.COM (Jim Adcock) (01/24/89)

/*******

For casual FFT'ers, here's a simple but elegant and efficient FFT
routine based on Rabiner and Gold pg 367

*******/

#include <complex.h>

// in-place fft a double precision complex array x of size (2 to the m)

void fft(complex* x, unsigned m)
{
    complex u, w, t;
    unsigned n = 1 << m;
    unsigned nv2 = n >> 1;
    unsigned nm1 = n - 1;
    unsigned j = 0;

    for (unsigned i = 0; i < nm1; ++i)
    {
        if (i < j) {t = x[j]; x[j] = x[i]; x[i] = t;}
        unsigned k = nv2;
        while ((k-1) < j) {j -= k; k >> 1;}
        j += k;
    }

    for (unsigned l = 1; l <= m; ++l)
    {
        unsigned le = 1 << l;
        unsigned le1 = le >> 1;
        complex u(1.0, 0.0);
        complex w(cos(M_PI/double(le1)),sin(M_PI/double(le1)));
        for (j = 0; j < le1; ++j)
        {
            for (i = j; i < n; i += le)
            {
                unsigned ip = i + le1;
                t = x[ip] * u;
                x[ip] = x[i] - t;
                x[i] += t;
            }
            u *= w;
        }
    }
}