This commit was generated by cvs2svn to compensate for changes in r58,

which included commits to RCS files with non-trunk default branches.

svn path=/trunk/; revision=59

1 files changed, 1001 insertions, 0 deletions

diff --git a/pd/src/d_fftroutine.c b/pd/src/d_fftroutine.c
new file mode 100644
index 00000000..7f0d2523
--- /dev/null
+++ b/pd/src/d_fftroutine.c
@@ -0,0 +1,1001 @@
+/*****************************************************************************/
+/*                                                                           */
+/* Fast Fourier Transform                                                    */
+/* Network Abstraction, Definitions                                          */
+/* Kevin Peterson, MIT Media Lab, EMS                                        */
+/* UROP - Fall '86                                                           */
+/* REV: 6/12/87(KHP) - To incorporate link list of different sized networks  */
+/*                                                                           */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* added debug option 5/91 brown@nadia                                       */
+/* change sign at AAA                                                        */
+/*                                                                           */
+/* Fast Fourier Transform                                                    */
+/* FFT Network Interaction and Support Modules                               */
+/* Kevin Peterson, MIT Media Lab, EMS                                        */
+/* UROP - Fall '86                                                           */
+/* REV: 6/12/87(KHP) - Generalized to one procedure call with typed I/O      */
+/*                                                                           */
+/*****************************************************************************/
+
+/* Overview:
+	
+   My realization of the FFT involves a representation of a network of
+   "butterfly" elements that takes a set of 'N' sound samples as input and
+   computes the discrete Fourier transform.  This network consists of a 
+   series of stages (log2 N), each stage consisting of N/2 parallel butterfly
+   elements. Consecutive stages are connected by specific, predetermined flow 
+   paths, (see Oppenheim, Schafer for details) and each butterfly element has
+   an associated multiplicative coefficient.
+
+   FFT NETWORK:
+   -----------	
+      ____    _    ____    _    ____    _    ____    _    ____
+  o--|    |o-| |-o|    |o-| |-o|    |o-| |-o|    |o-| |-o|    |--o
+     |reg1|  | |  |W^r1|  | |  |reg1|  | |  |W^r1|  | |  |reg1|
+     |    |  | |  |    |  | |  |    |  | |  |    |  | |  |    | .....
+     |    |  | |  |    |  | |  |    |  | |  |    |  | |  |    |  
+  o--|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|--o
+             | |          | |          | |          | |
+             | |          | |          | |          | |
+      ____   | |   ____   | |   ____   | |   ____   | |   ____ 
+  o--|    |o-| |-o|    |o-| |-o|    |o-| |-o|    |o-| |-o|    |--o
+     |reg2|  | |  |W^r2|  | |  |reg2|  | |  |W^r2|  | |  |reg2|
+     |    |  | |  |    |  | |  |    |  | |  |    |  | |  |    | .....
+     |    |  | |  |    |  | |  |    |  | |  |    |  | |  |    |
+  o--|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|o-| |-o|____|--o
+             | |          | |          | |          | |
+             | |          | |          | |          | |
+       :      :     :      :     :      :     :      :     :
+       :      :     :      :     :      :     :      :     :
+       :      :     :      :     :      :     :      :     :
+       :      :     :      :     :      :     :      :     :
+       :      :     :      :     :      :     :      :     :
+
+      ____   | |   ____   | |   ____   | |   ____   | |   ____ 
+  o--|    |o-| |-o|    |o-| |-o|    |o-| |-o|    |o-| |-o|    |--o
+     |reg |  | |  |W^r |  | |  |reg |  | |  |W^r |  | |  |reg |
+     | N/2|  | |  | N/2|  | |  | N/2|  | |  | N/2|  | |  | N/2| .....
+     |    |  | |  |    |  | |  |    |  | |  |    |  | |  |    |
+  o--|____|o-|_|-o|____|o-|_|-o|____|o-|_|-o|____|o-|_|-o|____|--o
+
+              ^            ^            ^            ^
+    Initial   |  Bttrfly   |   Rd/Wrt   |   Bttrfly  |   Rd/Wrt
+    Buffer    |            |  Register  |            |  Register
+              |____________|____________|____________|
+                                 |
+                                 |
+                            Interconnect
+			       Paths
+
+   The use of "in-place" computation permits one to use only one set of 
+   registers realized by an array of complex number structures.  To describe
+   the coefficients for each butterfly I am using a two dimensional array
+   (stage, butterfly) of complex numbers.  The predetermined stage connections
+   will be described in a two dimensional array of indicies.  These indicies 
+   will be used to determine the order of reading at each stage of the    
+   computation.  
+*/
+
+
+/*****************************************************************************/
+/* INCLUDE FILES                                                             */
+/*****************************************************************************/
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+    /* the following is needed only to declare pd_fft() as exportable in NT */
+#include "m_pd.h"
+
+/* some basic definitions */
+#ifndef BOOL
+#define        BOOL                    int
+#define        TRUE                    1
+#define        FALSE                   0
+#endif
+
+#define        SAMPLE float	/* data type used in calculation */
+
+#define        SHORT_SIZE              sizeof(short)
+#define        INT_SIZE                sizeof(int)
+#define        FLOAT_SIZE              sizeof(float)
+#define        SAMPLE_SIZE             sizeof(SAMPLE)
+#define        PNTR_SIZE               sizeof(char *)
+
+#define        PI                      3.1415927
+#define        TWO_PI                  6.2831854
+
+/* type definitions for I/O buffers */
+#define        REAL                    0          /* real only          */
+#define        IMAG                    2          /* imaginary only     */
+#define        RECT                    8          /* real and imaginary */
+#define        MAG                     16         /* magnitude only     */
+#define        PHASE                   32         /* phase only         */
+#define        POLAR                   64         /* magnitude and phase*/
+
+/* scale definitions for I/O buffers */
+#define        LINEAR                  0
+#define        DB                      1          /* 20log10            */
+
+/* transform direction definition */
+#define        FORWARD                 1          /* Forward FFT        */
+#define        INVERSE                 2          /* Inverse FFT        */
+
+/* window type definitions */
+#define        HANNING                 1
+#define        RECTANGULAR             0
+
+
+
+/* network structure definition */
+
+typedef struct Tfft_net {
+	int             n;
+	int             stages;
+	int             bps;
+	int             direction;
+	int             window_type;
+	int             *load_index;
+	SAMPLE          *window, *inv_window;
+        SAMPLE          *regr;
+	SAMPLE          *regi;
+	SAMPLE          **indexpr;
+	SAMPLE          **indexpi;
+	SAMPLE          **indexqr;
+	SAMPLE          **indexqi;
+        SAMPLE          *coeffr, *inv_coeffr;
+        SAMPLE          *coeffi, *inv_coeffi;
+	struct Tfft_net *next;  
+} FFT_NET;
+
+
+void cfft(int trnsfrm_dir, int npnt, int window,
+    float *source_buf, int source_form, int source_scale,
+    float *result_buf, int result_form, int result_scale, int debug);
+
+
+/*****************************************************************************/
+/* GLOBAL DECLARATIONS                                                       */
+/*****************************************************************************/
+
+static FFT_NET  *firstnet;
+
+/* prototypes */
+
+void net_alloc(FFT_NET *fft_net);
+void net_dealloc(FFT_NET *fft_net);
+int power_of_two(int n);
+void create_hanning(SAMPLE *window, int n, SAMPLE scale);
+void create_rectangular(SAMPLE *window, int n, SAMPLE scale);
+void short_to_float(short *short_buf, float *float_buf, int n);
+void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
+    int buf_scale, int trnsfrm_dir);
+void compute_fft(FFT_NET  *fft_net);
+void store_registers(FFT_NET	*fft_net, float *buf, int buf_form,
+    int buf_scale, int debug);
+void build_fft_network(FFT_NET *fft_net, int n, int window_type);
+
+/*****************************************************************************/
+/* GENERALIZED FAST FOURIER TRANSFORM MODULE                                 */
+/*****************************************************************************/
+
+void cfft(int trnsfrm_dir, int npnt, int window,
+    float *source_buf, int source_form, int source_scale,
+    float *result_buf, int result_form, int result_scale, int debug)
+
+/* modifies: result_buf
+   effects:  Computes npnt FFT specified by form, scale, and dir parameters.  
+         Source samples (single precision float) are taken from soure_buf and 
+	 the transfrmd representation is stored in result_buf (single precision
+	 float).  The parameters are defined as follows:
+	
+	 trnsfrm_dir = FORWARD | INVERSE
+	 npnt        = 2^k for some any positive integer k
+	 window      = HANNING | RECTANGULAR
+	 (RECT = real and imag parts, POLAR = magnitude and phase)
+	 source_form = REAL | IMAG | RECT | POLAR  
+	 result_form = REAL | IMAG | RECT | MAG | PHASE | POLAR
+	 xxxxxx_scale= LINEAR | DB ( 20log10 |mag| )
+	 
+	 The input/output buffers are stored in a form appropriate to the type.
+	 For example: REAL  => {real, real, real ...}, 
+	              MAG   => {mag, mag, mag, ... },
+		      RECT  => {real, imag, real, imag, ... },
+		      POLAR => {mag, phase, mag, phase, ... }.
+
+         To look at the magnitude (in db) of a 1024 point FFT of a real time 
+	 signal we have:
+
+	 fft(FORWARD, 1024, RECTANGULAR, input, REAL, LINEAR, output, MAG, DB)
+
+	 All possible input and output combinations are possible given the 
+	 choice of type and scale parameters.
+*/
+
+{
+         FFT_NET	 *thisnet = (FFT_NET *)0;
+	 FFT_NET         *lastnet = (FFT_NET *)0;
+	 
+	 /* A linked list of fft networks of different sizes is maintained to
+	    avoid building with every call.  The network is built on the first
+	    call but reused for subsequent calls requesting the same size 
+	    transformation.
+	 */
+   
+	 thisnet=firstnet;
+	 while (thisnet) {
+	     if (!(thisnet->n == npnt) || !(thisnet->window_type == window)) { 
+	       /* current net doesn't match size or window type */
+	       lastnet=thisnet;
+	       thisnet=thisnet->next;
+	       continue;                  /* keep looking */
+	     }
+
+	     else {                       /* network matches desired size */
+	       load_registers(thisnet, source_buf, source_form, source_scale, 
+			      trnsfrm_dir);
+	       compute_fft(thisnet);      /* do transformation */
+	       store_registers(thisnet, result_buf, result_form, result_scale,debug);
+	       return;
+	     }
+         }	     
+
+	 /* none of existing networks match required size*/
+
+	 if (lastnet) {                 /* add new network to end of list */
+	   thisnet = (FFT_NET *)malloc(sizeof(FFT_NET));      /* allocate */
+	   thisnet->next = 0;
+	   lastnet->next = thisnet;     /* add to end of list             */
+	 }
+	 else {                         /* first network to be created    */
+	   thisnet=firstnet=(FFT_NET *)malloc(sizeof(FFT_NET)); /* alloc. */
+	   thisnet->next = 0;
+	 }
+
+	 /* build new network and compute transformation */
+	 build_fft_network(thisnet, npnt, window);
+	 load_registers(thisnet, source_buf, source_form, source_scale, 
+			trnsfrm_dir);
+	 compute_fft(thisnet);
+	 store_registers(thisnet, result_buf, result_form, result_scale,debug);
+	 return;
+}
+
+void fft_clear(void)
+
+/* effects: Deallocates all preserved FFT networks.  Should be used when 
+         finished with all computations.
+*/
+
+{
+         FFT_NET      *thisnet, *nextnet;
+
+	 if (firstnet) {
+	   thisnet=firstnet;
+	   do {
+	     nextnet = thisnet->next;
+	     net_dealloc(thisnet);
+	     free((char *)thisnet);
+	   } while (thisnet = nextnet);
+	 }
+}
+	   
+
+/*****************************************************************************/
+/* NETWORK CONSTRUCTION                                                      */
+/*****************************************************************************/
+
+void build_fft_network(FFT_NET *fft_net, int n, int window_type)
+
+
+/* modifies:fft_net
+   effects: Constructs the fft network as described in fft.h.  Butterfly
+         coefficients, read/write indicies, bit reversed load indicies,
+	 and array allocations are computed.
+*/
+
+{
+	 int cntr, i, j, s; 
+	 int       stages, bps;
+         int       **p, **q, *pp, *qp;
+	 SAMPLE     two_pi_div_n = TWO_PI / n;
+
+
+	 /* network definition */
+	 fft_net->n   = n;
+	 fft_net->bps = bps = n/2;
+	 for (i = 0, j = n; j > 1; j >>= 1, i++);
+	 fft_net->stages = stages = i;
+	 fft_net->direction = FORWARD;
+	 fft_net->window_type = window_type;
+	 fft_net->next = (FFT_NET *)0;
+
+	 /* allocate registers, index, coefficient arrays */
+	 net_alloc(fft_net);
+
+
+	 /* create appropriate windows */
+	 if (window_type==HANNING)   {
+		  create_hanning(fft_net->window, n, 1.);
+		  create_hanning(fft_net->inv_window, n, 1./n);
+	 }
+	 else {
+		  create_rectangular(fft_net->window, n, 1.);
+		  create_rectangular(fft_net->inv_window, n, 1./n);
+	 }
+
+
+	 /* calculate butterfly coefficients */ {
+		  
+		  int       num_diff_coeffs, power_inc, power;
+		  SAMPLE *coeffpr     = fft_net->coeffr;
+		  SAMPLE *coeffpi     = fft_net->coeffi;
+		  SAMPLE *inv_coeffpr = fft_net->inv_coeffr;
+		  SAMPLE *inv_coeffpi = fft_net->inv_coeffi;
+		  
+		  /* stage one coeffs are 1 + 0j */
+		  for (i = 0; i < bps; i++) {
+			   *coeffpr = *inv_coeffpr = 1.;
+			   *coeffpi = *inv_coeffpi = 0.;
+			   coeffpr++; inv_coeffpr++;
+			   coeffpi++; inv_coeffpi++;
+		  }
+
+		  /* stage 2 to last stage coeffs need calculation */
+		  /* (1<<r <=> 2^r */
+		  for (s = 2; s <= stages; s++) {
+			   
+			   num_diff_coeffs = n / (1 << (stages - s + 1)); 
+			   power_inc       = 1 << (stages -s);
+			   cntr            = 0;
+
+			   for (i = bps/num_diff_coeffs; i > 0; i--) {
+
+		              power  = 0;
+
+			      for (j = num_diff_coeffs; j > 0; j--) {
+			         *coeffpr     = cos(two_pi_div_n*power);
+				 *inv_coeffpr = cos(two_pi_div_n*power);
+/* AAA change these signs */	 *coeffpi     = -sin(two_pi_div_n*power);
+/* change back */		 *inv_coeffpi = sin(two_pi_div_n*power);
+				 power += power_inc;
+				 coeffpr++; inv_coeffpr++;
+				 coeffpi++; inv_coeffpi++;
+			      }
+		           }
+		  }
+	 }
+
+	 /* calculate network indicies:  stage exchange indicies are 
+	    calculated and then used as offset values from the base
+	    register locations.  The final addresses are then stored in
+	    fft_net.
+	 */ {
+
+		  int       index, inc;
+		  SAMPLE **indexpr = fft_net->indexpr;
+		  SAMPLE **indexpi = fft_net->indexpi;
+		  SAMPLE **indexqr = fft_net->indexqr;
+		  SAMPLE **indexqi = fft_net->indexqi;
+		  SAMPLE *regr     = fft_net->regr;
+		  SAMPLE *regi     = fft_net->regi; 
+
+
+		  /* allocate temporary 2d stage exchange index, 1d temp 
+		     load index */
+		  p = (int **)malloc(stages * PNTR_SIZE);
+		  q = (int **)malloc(stages * PNTR_SIZE);
+
+		  for (s = 0; s < stages; s++) {
+			   p[s] = (int *)malloc(bps * INT_SIZE);
+			   q[s] = (int *)malloc(bps * INT_SIZE);
+		  }
+
+		  /* calculate stage exchange indicies: */
+		  for (s = 0; s < stages; s++) {
+			   pp = p[s];
+			   qp = q[s];
+			   inc    = 1 << s;
+			   cntr   = 1 << (stages-s-1);
+			   i      = j = index = 0;
+
+			   do {
+				    do {
+					     qp[i]   = index + inc;
+					     pp[i++] = index++;
+				    }  while (++j < inc);
+				    index = qp[i-1] + 1;
+				    j = 0;
+			   }        while (--cntr);
+		  }
+
+		  /* compute actual address values using indicies as offsets */
+		  for (s = 0; s < stages; s++) {
+			   for (i = 0; i < bps; i++) {
+				    *indexpr++ = regr + p[s][i];
+				    *indexpi++ = regi + p[s][i];
+				    *indexqr++ = regr + q[s][i];
+				    *indexqi++ = regi + q[s][i];
+			   }
+		  }
+	 }
+
+
+         /* calculate load indicies (bit reverse ordering) */
+         /* bit reverse ordering achieved by passing normal
+	    order indicies backwards through the network */
+		  
+	 /* init to normal order indicies */ {
+		  int *load_index,*load_indexp;
+		  int *temp_indexp, *temp_index;
+		  temp_index=temp_indexp=(int *)malloc(n * INT_SIZE);
+			   
+		  i = 0; j = n;
+		  load_index = load_indexp = fft_net->load_index;
+			   
+		  while (j--)
+		    *load_indexp++ = i++;
+
+	/* pass indicies backwards through net */
+		  for (s = stages - 1; s > 0; s--) {
+			   pp = p[s];
+			   qp = q[s];
+
+			   for (i = 0; i < bps; i++) {
+				    temp_index[pp[i]]=load_index[2*i];
+				    temp_index[qp[i]]=load_index[2*i+1];
+			   }
+			   j = n;
+			   load_indexp = load_index;
+			   temp_indexp = temp_index;
+			   while (j--) 
+			     *load_indexp++ = *temp_indexp++;
+		  }
+		  
+		  /* free all temporary arrays */
+		  free((char *)temp_index);
+		  for (s = 0; s < stages; s++) {
+			   free((char *)p[s]);free((char *)q[s]);
+		  }
+		  free((char *)p);free((char *)q);
+	 }
+}
+
+
+
+/*****************************************************************************/
+/* REGISTER LOAD AND STORE                                                   */
+/*****************************************************************************/
+
+void load_registers(FFT_NET *fft_net, float *buf, int buf_form,
+    int buf_scale, int trnsfrm_dir)
+
+/* effects:  Multiplies the input buffer with the appropriate window and
+         stores the resulting values in the initial registers of the
+	 network.  Input buffer must contain values appropriate to form.  
+	 For RECT, the buffer contains real num. followed by imag num, 
+	 and for POLAR, it contains magnitude followed by phase.  Pure
+	 inputs are listed normally.  Both LINEAR and DB scales are 
+	 interpreted.
+*/
+
+{
+	 int      *load_index = fft_net->load_index;
+	 SAMPLE *window;
+	 int index, i = 0, n = fft_net->n;
+
+	 if      (trnsfrm_dir==FORWARD)   window = fft_net->window;
+	 else if (trnsfrm_dir==INVERSE)   window = fft_net->inv_window;
+	 else {
+		  fprintf(stderr, "load_registers:illegal transform direction\n"); 
+		  exit(0);
+	 }
+	 fft_net->direction = trnsfrm_dir;
+
+	 switch(buf_scale) {
+	 case LINEAR: {
+
+	   switch (buf_form) {
+	   case REAL: {                    /* pure REAL */
+	     while (i < fft_net->n) {  
+	       index = load_index[i];
+	       fft_net->regr[i]=(SAMPLE)buf[index] * window[index];
+	       fft_net->regi[i]=0.;
+	       i++;			                       
+	     }
+	   } break;
+
+	   case IMAG: {                    /* pure IMAGinary */
+	     while (i < fft_net->n) {  
+	       index = load_index[i];
+	       fft_net->regr[i]=0;
+	       fft_net->regi[i]=(SAMPLE)buf[index] * window[index];
+	       i++;			       
+	     }                
+	   } break;
+
+           case RECT: {                    /* both REAL and IMAGinary */
+	     while (i < fft_net->n) {
+	       index = load_index[i];
+	       fft_net->regr[i]=(SAMPLE)buf[index*2]   * window[index];
+	       fft_net->regi[i]=(SAMPLE)buf[index*2+1] * window[index];
+	       i++;
+	     }
+	   } break;	 
+  	 
+	   case POLAR: {                   /* magnitude followed by phase */
+	     while (i < fft_net->n) {
+	       index = load_index[i];
+	       fft_net->regr[i]=(SAMPLE)(buf[index*2] * cos(buf[index*2+1])) 
+		                                      * window[index];
+	       fft_net->regi[i]=(SAMPLE)(buf[index*2] * sin(buf[index*2+1])) 
+	                                              * window[index];
+	       i++;			       
+	     }                
+	   } break;
+
+	   default: {
+	     fprintf(stderr, "load_registers:illegal input form\n"); 
+	     exit(0);
+	   } break;
+           }
+	 } break;
+
+         case DB: {
+	  
+	   switch (buf_form) {
+	   case REAL: {                     /* log pure REAL */
+	     while (i < fft_net->n) {  
+	       index = load_index[i];
+	       fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index]) 
+		 * window[index];    /* window scaling after linearization */
+	       fft_net->regi[i]=0.;
+	       i++;			                       
+	     }
+	   } break;
+
+	   case IMAG: {                     /* log pure IMAGinary */
+	     while (i < fft_net->n) {  
+	       index = load_index[i];
+	       fft_net->regr[i]=0.;
+	       fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index])
+		    * window[index];
+	       i++;			       
+	    }                
+	   } break;
+
+           case RECT: {                     /* log REAL and log IMAGinary */
+	     while (i < fft_net->n) {
+	       index = load_index[i];
+	       fft_net->regr[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2])
+		 * window[index];
+	       fft_net->regi[i]=(SAMPLE)pow(10., (1./20.)*buf[index*2+1]) 
+		 * window[index];
+	       i++;
+	     }
+	   } break;	 
+  	 
+	   case POLAR: {                    /* log mag followed by phase */
+	     while (i < fft_net->n) {
+	       index = load_index[i];
+	       fft_net->regr[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
+			     * cos(buf[index*2+1])) * window[index];
+	       fft_net->regi[i]=(SAMPLE)(pow(10., (1./20.)*buf[index*2])
+			     * sin(buf[index*2+1])) * window[index];
+	       i++;			       
+	     }                
+	   } break;
+
+	   default: {
+	     fprintf(stderr, "load_registers:illegal input form\n"); 
+	     exit(0);
+	   } break;
+	   }
+	 } break;
+
+         default: {
+	   fprintf(stderr, "load_registers:illegal input scale\n"); 
+	   exit(0);
+	 } break;
+	 }
+}
+
+
+void store_registers(FFT_NET	*fft_net, float *buf, int buf_form,
+    int buf_scale, int debug)
+
+/* modifies: buf
+   effects:  Writes the final contents of the network registers into buf in 
+         either linear or db scale, polar or rectangular form.  If any of 
+	 the pure forms(REAL, IMAG, MAG, or PHASE) are used then only the 
+	 corresponding part of the registers is stored in buf.
+*/
+
+{
+	 int	    i;
+	 SAMPLE	    real, imag, mag, phase;
+	 int	    n;
+
+	 i = 0;
+	 n = fft_net->n;
+
+	 switch (buf_scale) {
+	 case LINEAR: {
+
+	   switch (buf_form) {
+	   case REAL: {                        /* pure REAL */
+	     do {
+	       *buf++ = (float)fft_net->regr[i];
+	     } while (++i < n);  
+	   } break;
+
+	   case IMAG: {                        /* pure IMAGinary */
+	     do {
+	       *buf++ = (float)fft_net->regi[i];
+	     } while (++i < n);  
+	   } break;
+
+           case RECT: {	         	       /* both REAL and IMAGinary */   
+	     do {
+	       *buf++ = (float)fft_net->regr[i];
+	       *buf++ = (float)fft_net->regi[i];
+	     } while (++i < n);  
+	   } break;
+
+	   case MAG: {                         /* magnitude only */
+	     do {
+	       real  = fft_net->regr[i];
+	       imag  = fft_net->regi[i];
+	       *buf++ = (float)sqrt(real*real+imag*imag);
+	     } while (++i < n);
+	   } break;
+
+	   case PHASE: {                       /* phase only */
+	     do {
+	       real  = fft_net->regr[i];
+	       imag  = fft_net->regi[i];
+	       if (real > .00001) 
+		 *buf++ = (float)atan2(imag, real);
+	       else {                          /* deal with bad case */
+		 if (imag > 0){      *buf++ = PI / 2.;
+		     if(debug) fprintf(stderr,"real=0 and imag > 0\n");}
+		 else if (imag < 0){ *buf++ = -PI / 2.;
+		     if(debug) fprintf(stderr,"real=0 and imag < 0\n");}
+		 else {              *buf++ = 0;
+		     if(debug) fprintf(stderr,"real=0 and imag=0\n");}
+	       }
+	     } while (++i < n);
+	   } break;
+
+	   case POLAR: {                       /* magnitude and phase */
+	     do {
+	       real    = fft_net->regr[i];
+	       imag    = fft_net->regi[i];
+	       *buf++  = (float)sqrt(real*real+imag*imag);
+	       if (real)                       /* a hack to avoid div by zero */
+		 *buf++ = (float)atan2(imag, real);
+	       else {                          /* deal with bad case */
+		 if (imag > 0)      *buf++ = PI / 2.;
+		 else if (imag < 0) *buf++ = -PI / 2.;
+		 else               *buf++ = 0;
+	       }
+	     } while (++i < n);
+	   } break;
+
+	   default: {
+	     fprintf(stderr, "store_registers:illegal output form\n");
+	     exit(0);
+	   } break;
+	   }
+	 } break;
+		
+         case DB: {
+
+	   switch (buf_form) {
+	   case REAL: {                        /* real only */
+	     do {
+	       *buf++ = (float)20.*log10(fft_net->regr[i]);
+	     } while (++i < n);
+	   } break;
+
+	   case IMAG: {                        /* imag only */
+	     do {
+	       *buf++ = (float)20.*log10(fft_net->regi[i]);
+	     } while (++i < n);
+	   } break;
+
+	   case RECT: {                        /* real and imag */
+	     do {
+	       *buf++ = (float)20.*log10(fft_net->regr[i]);
+	       *buf++ = (float)20.*log10(fft_net->regi[i]);
+	     } while (++i < n);  
+	   } break;
+
+	   case MAG: {                         /* magnitude only  */
+	     do {
+	       real  = fft_net->regr[i];
+	       imag  = fft_net->regi[i];
+	       *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));  
+	     } while (++i < n);
+	   } break;
+
+	   case PHASE: {                       /* phase only */
+	     do {
+	       real  = fft_net->regr[i];
+	       imag  = fft_net->regi[i];
+	       if (real) 
+		 *buf++ = (float)atan2(imag, real);
+	       else {                          /* deal with bad case */
+		 if (imag > 0)      *buf++ = PI / 2.;
+		 else if (imag < 0) *buf++ = -PI / 2.;
+		 else               *buf++ = 0;
+	       }
+	     } while (++i < n);
+	   } break;
+
+	   case POLAR: {                       /* magnitude and phase */
+	     do {
+	       real  = fft_net->regr[i];
+	       imag  = fft_net->regi[i];
+	       *buf++ = (float)20.*log10(sqrt(real*real+imag*imag));  	       
+	       if (real) 
+		 *buf++ = (float)atan2(imag, real);
+	       else {                          /* deal with bad case */
+		 if (imag > 0)      *buf++ = PI / 2.;
+		 else if (imag < 0) *buf++ = -PI / 2.;
+		 else               *buf++ = 0;
+	       }
+	     } while (++i < n);
+	   } break;
+
+  	   default: {
+	     fprintf(stderr, "store_registers:illegal output form\n");
+	     exit(0);
+	   } break;
+	   } 
+	 } break;
+
+         default: {
+	   fprintf(stderr, "store_registers:illegal output scale\n");
+	   exit(0);
+	 } break;
+         }
+}
+
+
+
+/*****************************************************************************/
+/* COMPUTE TRANSFORMATION                                                    */
+/*****************************************************************************/
+
+void compute_fft(FFT_NET  *fft_net)
+
+
+/* modifies: fft_net
+   effects: Passes the values (already loaded) in the registers through
+	 the network, multiplying with appropriate coefficients at each 
+	 stage.  The fft result will be in the registers at the end of
+	 the computation.  The direction of the transformation is indicated
+	 by the network flag 'direction'.  The form of the computation is:
+
+	 X(pn) = X(p) + C*X(q)
+	 X(qn) = X(p) - C*X(q)
+
+	 where X(pn,qn) represents the output of the registers at each stage.  
+	 The calculations are actually done in place.  Register pointers are 
+	 used to speed up the calculations.
+
+	 Register and coefficient addresses involved in the calculations 
+	 are stored sequentially and are accessed as such. fft_net->indexp,
+	 indexq contain pointers to the relevant addresses, and fft_net->coeffs, 
+	 inv_coeffs points to the appropriate coefficients at each stage of the 
+	 computation.
+*/
+
+{
+	 SAMPLE	    **xpr, **xpi, **xqr, **xqi, *cr, *ci;
+	 int	    i;
+	 SAMPLE	    tpr, tpi, tqr, tqi;
+	 int	    bps = fft_net->bps;
+	 int	    cnt = bps * (fft_net->stages - 1);
+
+	 /* predetermined register addresses and coefficients */
+	 xpr = fft_net->indexpr;              
+	 xpi = fft_net->indexpi;              
+	 xqr = fft_net->indexqr;
+	 xqi = fft_net->indexqi;
+
+	 if (fft_net->direction==FORWARD) { 	/* FORWARD FFT coefficients */
+		  cr  = fft_net->coeffr;
+		  ci  = fft_net->coeffi;
+	 }
+	 else {     	    	    	    	/* INVERSE FFT coefficients */
+		  cr = fft_net->inv_coeffr;
+		  ci = fft_net->inv_coeffi;
+	 }
+
+	 /* stage one coefficients are 1 + 0j so C*X(q)=X(q)  */
+	 /* bps mults can be avoided                          */
+
+	 for (i = 0; i < bps; i++) {
+
+		  /* add X(p) and X(q) */
+		  tpr = **xpr + **xqr;
+		  tpi = **xpi + **xqi;
+		  tqr = **xpr - **xqr;
+		  tqi = **xpi - **xqi;
+		  
+		  /* exchange register with temp */
+		  **xpr = tpr;
+		  **xpi = tpi;
+		  **xqr = tqr;
+		  **xqi = tqi;
+
+		  /* next set of register for calculations: */
+		  xpr++; xpi++; xqr++; xqi++; cr++; ci++;
+
+	 }
+
+	 for (i = 0; i < cnt; i++) {
+		  
+		  /* mult X(q) by coeff C */
+		  tqr = **xqr * *cr - **xqi * *ci;
+		  tqi = **xqr * *ci + **xqi * *cr;
+
+		  /* exchange register with temp */
+		  **xqr = tqr;
+		  **xqi = tqi;
+
+		  /* add X(p) and X(q) */
+		  tpr = **xpr + **xqr;
+		  tpi = **xpi + **xqi;
+		  tqr = **xpr - **xqr;
+		  tqi = **xpi - **xqi;
+		  
+		  /* exchange register with temp */
+		  **xpr = tpr;
+		  **xpi = tpi;
+		  **xqr = tqr;
+		  **xqi = tqi;
+		  /* next set of register for calculations: */
+		  xpr++; xpi++; xqr++; xqi++; cr++; ci++;
+	 }
+}
+
+
+/****************************************************************************/
+/* SUPPORT MODULES                                                          */
+/****************************************************************************/
+
+void net_alloc(FFT_NET *fft_net)
+
+
+/* effects: Allocates appropriate two dimensional arrays and assigns
+	   correct internal pointers.
+*/
+
+{
+
+	 int      stages, bps, n;
+
+	 n      = fft_net->n;
+	 stages = fft_net->stages;
+	 bps    = fft_net->bps;
+
+
+ 	 /* two dimensional arrays with elements stored sequentially */
+
+	 fft_net->load_index  = (int *)malloc(n * INT_SIZE);
+	 fft_net->regr        = (SAMPLE *)malloc(n * SAMPLE_SIZE);
+	 fft_net->regi        = (SAMPLE *)malloc(n * SAMPLE_SIZE);
+	 fft_net->coeffr      = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
+	 fft_net->coeffi      = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
+	 fft_net->inv_coeffr  = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
+	 fft_net->inv_coeffi  = (SAMPLE *)malloc(stages*bps*SAMPLE_SIZE);
+	 fft_net->indexpr     = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
+	 fft_net->indexpi     = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
+	 fft_net->indexqr     = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
+	 fft_net->indexqi     = (SAMPLE **)malloc(stages * bps * PNTR_SIZE);
+
+	 /* one dimensional load window */
+	 fft_net->window      = (SAMPLE *)malloc(n * SAMPLE_SIZE);
+	 fft_net->inv_window  = (SAMPLE *)malloc(n * SAMPLE_SIZE);
+}
+
+void net_dealloc(FFT_NET *fft_net)
+
+
+/* effects: Deallocates given FFT network.
+*/
+
+{
+
+	 free((char *)fft_net->load_index);  
+	 free((char *)fft_net->regr);        
+	 free((char *)fft_net->regi);        
+	 free((char *)fft_net->coeffr);      
+	 free((char *)fft_net->coeffi);      
+	 free((char *)fft_net->inv_coeffr);  
+	 free((char *)fft_net->inv_coeffi);  
+	 free((char *)fft_net->indexpr);     
+	 free((char *)fft_net->indexpi);     
+	 free((char *)fft_net->indexqr);   
+	 free((char *)fft_net->indexqi);   
+	 free((char *)fft_net->window);
+	 free((char *)fft_net->inv_window);
+}
+
+
+BOOL power_of_two(n)
+
+int               n;
+
+/* effects: Returns TRUE if n is a power of two, otherwise FALSE.
+*/
+
+{
+	 int      i;
+
+	 for (i = n; i > 1; i >>= 1) 
+		  if (i & 1) return FALSE;        /* more than one bit high */
+	 return TRUE;
+}
+
+
+void create_hanning(SAMPLE *window, int n, SAMPLE scale)
+
+/* effects: Fills the buffer window with a hanning window of the appropriate
+         size scaled by scale.
+*/
+
+{
+	 SAMPLE	    a, pi_div_n = PI/n;
+	 int	    k;
+
+	 for (k=1; k <= n; k++) {
+		  a = sin(k * pi_div_n);
+		  *window++ = scale * a * a;
+	 }
+}
+
+
+void create_rectangular(SAMPLE *window, int n, SAMPLE scale)
+
+/* effects: Fills the buffer window with a rectangular window of the
+   appropriate size of height scale.
+*/
+
+{
+	 while (n--)
+	   *window++ = scale;
+}
+
+
+void short_to_float(short *short_buf, float *float_buf, int n)
+
+/* effects; Converts short_buf to floats and stores them in float_buf.
+*/
+
+{
+	 while (n--) {
+		  *float_buf++ = (float)*short_buf++;
+	 }
+}
+
+
+/* here's the meat: */
+
+void pd_fft(float *buf, int npoints, int inverse)
+{
+  double renorm;
+  float *fp, *fp2;
+  int i;
+  renorm = (inverse ? npoints : 1.);
+  cfft((inverse ? INVERSE : FORWARD), npoints, RECTANGULAR, 
+       buf, RECT, LINEAR, buf, RECT, LINEAR, 0);
+  for (i = npoints << 1, fp = buf; i--; fp++) *fp *= renorm;
+}