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