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 }