00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
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
00034
00035
00036
00037 void fftPrepDebug( int level )
00038 {
00039 Debug = level;
00040 return;
00041 }
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 long buildFacList( long new_nfft )
00057 {
00058 long i, j, next_nfft;
00059 FACT *this, *next, *ins_before;
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
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
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
00094 this = fft_list;
00095 while (this->nfft < last_nfft )
00096 {
00097
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
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
00127
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 {
00134 continue;
00135 }
00136
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
00144 next->next = ins_before->next;
00145 ins_before->next = next;
00146 }
00147
00148 this = this->next;
00149 last_nfft = this->nfft;
00150 }
00151 return 0;
00152 }
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169 long prepFFT( long n, FACT **pf )
00170 {
00171 FACT *this;
00172 int rc;
00173
00174
00175 if ( n > last_nfft)
00176 {
00177 if ( (rc = buildFacList( n )) < 0 )
00178 return rc;
00179 }
00180
00181
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
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
00234
00235
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 }