Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

fft_prep.c

Go to the documentation of this file.
00001 /*
00002  *   THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE
00003  *   CHECKED IT OUT USING THE COMMAND CHECKOUT.
00004  *
00005  *    $Id: 
00006  *
00007  *    Revision history:
00008  *     $Log:
00009  *
00010  *
00011  *
00012  */
00013 
00014 /*
00015  * fft_prep.c: routines for setting up the structures needed for the 
00016  * Temperton FFT99 package.
00017  */
00018 
00019 #include <stdlib.h>
00020 #include <stdio.h>
00021 #include <fft_prep.h>
00022 #include <fft99.h>
00023 #include <earthworm.h>
00024 
00025 static long last_nfft = 1;
00026 static FACT *fft_list = (FACT *)NULL;
00027 static int prime[3] = {2, 3, 5};
00028 static int Debug = 0;
00029 
00030 static int makeTrigs( FACT *this );
00031 
00032 /*
00033  * fftPrepDebug: set the debug level for all functions in fft_prep.c
00034  *        Debug levels: 0 no debug output
00035  *                      1 dump FAClist on logic errors.
00036  */
00037 void fftPrepDebug( int level )
00038 {
00039   Debug = level;
00040   return;
00041 }
00042 
00043 
00044 /* buildFacList: build a linked list of FFT factor structures. This list is
00045  * used to decide the size of FFT to use. The FFT routines can handle 
00046  * sizes that are multiples of powers of 2, 3, and 5. This provides greater
00047  * flexibility than sizes that are only powers of 2.
00048  *
00049  * argument: new_nfft: the largest fft size required for this list.
00050  *   This function may be called again to extend the list further.
00051  * returns: 0 on success
00052  *         -1 on failure (out of memory; not logged here)
00053  *         -2 on logic errors (logged here)
00054  */
00055 
00056 long buildFacList( long new_nfft )
00057 {
00058   long i, j, next_nfft;
00059   FACT *this, *next, *ins_before;
00060 
00061   /*
00062    * The idea is to start at `last_nfft' and multiply it by 2, 3, and 5, 
00063    * creating a new FACT element on the list for each one. After this is
00064    * done, then we advance on the list ONE element (not the three that we 
00065    * just created. From this next list element, we set `last_nfft' to the
00066    * nfft value of the element. In this way, the list contains an ordered
00067    * set of all the multiples of powers of 2, 3, and 5.
00068    */
00069 
00070   /* First time through is a special case: there are no elements on the 
00071    * list, so we have to create one. We'll assume that new_nfft > 5
00072    * though it really doesn't matter */
00073   if (last_nfft == 1)
00074   {
00075     if ( (fft_list = (FACT *)calloc(1, sizeof(FACT))) == (FACT *)NULL)
00076       return -1;
00077     this = fft_list;
00078     this->nfft = prime[0];
00079     this->fact_power[0] = 1;
00080     next = this;
00081     for (i = 1; i < N_RADIX; i++)
00082     {
00083       /* Get the next factor structure and fill it in */
00084       if ( (next->next = (FACT *)calloc(1, sizeof(FACT))) == (FACT *)NULL)
00085         return -1;
00086       next = next->next;
00087       next->nfft = last_nfft * prime[i];
00088       next->fact_power[i]++;
00089     }
00090   } 
00091   else
00092   {
00093     /* Walk the list to the last (newest, largest) entry */
00094     this = fft_list;
00095     while (this->nfft < last_nfft )
00096     {
00097       /* DEBUG; shouldn't happen if logic is correct */
00098       if (this->next == (FACT *) NULL)
00099       {
00100         logit("et", "buildFacList: Walked off end of list this %ld last %ld\n",
00101               this->nfft, last_nfft);
00102         if (Debug > 0)
00103           printFacList();
00104         return -2;
00105       }
00106       this = this->next;
00107     }
00108     if (this->nfft != last_nfft)
00109     {
00110       logit("et", "buildFacList: last %ld does not match this %ld\n", 
00111             last_nfft, this->nfft);
00112       if (Debug > 0)
00113         printFacList();
00114       return -2;
00115     }
00116   }
00117   last_nfft = this->nfft;
00118   
00119   /* Add new entries to the list up to the requested "new_nfft" */
00120   while (this->nfft <= new_nfft)
00121   {
00122     ins_before = this;
00123     for (i = 0; i < N_RADIX; i++)
00124     {
00125       next_nfft = last_nfft * prime[i];
00126       /* Walk the list to see where the next element will go, or if it 
00127        * is a dupe */
00128       while( ins_before->next != (FACT *)NULL && 
00129              ins_before->next->nfft < next_nfft)
00130         ins_before = ins_before->next;
00131       if (ins_before->next != (FACT *)NULL && 
00132           ins_before->next->nfft == next_nfft)
00133       {  /* Found a dupe; skip it */
00134         continue;
00135       }
00136       /* Get the next factor structure and fill it in */
00137       if ( ( next = (FACT *)calloc(1, sizeof(FACT))) == (FACT *)NULL)
00138         return -1;
00139       next->nfft = next_nfft;
00140       for (j = 0; j < N_RADIX; j++)
00141         next->fact_power[j] = this->fact_power[j];
00142       next->fact_power[i]++;
00143       /* Insert in the list */
00144       next->next = ins_before->next;
00145       ins_before->next = next;
00146     }
00147     /* don't need current factor anymore, move to next one */
00148     this = this->next;
00149     last_nfft = this->nfft;
00150   }
00151   return 0;
00152 }
00153 
00154 
00155 /*
00156  * prepFFT: find a `suitable' FFT size for proposed size n, and return
00157  * a pointer to the FACT structure that has been prepared for this FFT.
00158  * Currently `suitable' is defined as the next even nfft value larger than n.
00159  * In the future, a more intelligent sizing algorithm may be employed.
00160  *
00161  * Arguments: n   the proposed FFT size
00162  *           *pf  pointer to the FACT structure that will be filled in
00163  *                by prepFFT for this FFT.
00164  * returns:  size of FFT on success, 
00165  *           -1 on failure (out of memory; not logged)
00166  *           -2 on other errors, logged here
00167  */
00168 
00169 long prepFFT( long n, FACT **pf )
00170 {
00171   FACT *this;
00172   int rc;
00173   
00174   /* Make sure we have enough entries in the list of FACT structures */
00175   if ( n > last_nfft)
00176   {
00177     if ( (rc = buildFacList( n )) < 0 ) 
00178       return rc;
00179   }
00180   
00181   /* Walk the list, looking for `suitable', even nfft value */
00182   this = fft_list;
00183   while( (this->nfft < n || this->nfft % 2 == 1) && 
00184          this->next != (FACT *)NULL )
00185     this = this->next;
00186 
00187   if ( this->trigs == (double *)NULL)
00188   {
00189     if (makeTrigs( this ) < 0)
00190       return -1;
00191   }
00192   
00193   *pf = this;
00194   return this->nfft;
00195 }
00196 
00197   
00198 void printFacList( )
00199 {
00200   FACT *this;
00201 
00202   this = fft_list;
00203   
00204   while( this != (FACT *)NULL)
00205   {
00206     logit("", "%6ld %2ld %2ld %2ld\n", this->nfft, this->fact_power[0], 
00207            this->fact_power[1], this->fact_power[2]);
00208     this = this->next;
00209   }
00210   return;
00211 }
00212 
00213 /* trimFacList: remove unwanted entries from the beginning of the FFT list. */
00214 void trimFacList( long n )
00215 {
00216   FACT *this;
00217   
00218   if (n >= last_nfft)
00219     return;
00220   
00221   this = fft_list;
00222   
00223   while(this->nfft < n)
00224   {
00225     fft_list = this->next;
00226     free( this );
00227     this = fft_list;
00228   }
00229   return;
00230 }
00231 
00232 /*
00233  * makeTrigs: Set up the trigs and ifax arrays for fft99. 
00234  *    Returns: 0 on success
00235  *            -1 when out of memory
00236  */
00237 static int makeTrigs( FACT *this )
00238 {
00239   long ntrigs;
00240   
00241   ntrigs = 3 * this->nfft / 2 + 1;
00242   if ( (this->trigs = (double *)calloc(ntrigs, sizeof(double))) ==
00243        (double *)NULL)
00244     return -1;
00245   if ( (this->ifax = (long *)calloc(13, sizeof(long))) == (long *)NULL)
00246   {
00247     free(this->trigs);
00248     this->trigs = (double *)NULL;
00249     return -1;
00250   }
00251   fftfax(this->nfft, this->ifax, this->trigs);
00252 
00253   return 0;
00254 } 

Generated on Tue May 6 09:16:01 2003 for Earthworm Libs by doxygen1.3-rc3