142 {
143 char *indepQuantity, **depenQuantity, **exclude, **depenQuantityPair;
144 long depenQuantities, excludes;
145 char *input, *output;
146 long iArg, j, readCode, noWarnings, items;
147 int64_t i, rows, rowsToUse;
148 unsigned long flags, pairFlags, tmpFlags, pipeFlags, majorOrderFlag;
149 SCANNED_ARG *scArg;
151 double *tdata, *data, t0, dt;
152 double fracRMSChangeLimit, fracFreqAccuracyLimit;
153 int32_t frequenciesDesired, maxFrequencies, freqCycleLimit;
154 short truncate;
155 double *frequency, *amplitude = NULL, *phase = NULL, *significance = NULL, *phase1 = NULL, *amplitude1 = NULL, *significance1 = NULL;
156 long *frequencyIndex, *amplitudeIndex, *phaseIndex, *significanceIndex, pairs;
157 long *amplitudeIndex1, *phaseIndex1, *significanceIndex1;
158 short columnMajorOrder = -1;
159
161
162#ifdef DEBUG
163 if (1) {
164 long code;
165 double x, y;
166 x = 1.1;
167 code = OneDFunctionOptimize(&y, &x, 0.07, -4, 4, trialFn, 50, 1e-6, 0, 1);
168 fprintf(stderr, "code: %ld x=%e, y=%e\n", code, x, y);
169
170 x = .9;
171 code = OneDFunctionOptimize(&y, &x, 0.15, -4, 4, trialFn, 50, 1e-6, 0, 1);
172 fprintf(stderr, "code: %ld x=%e, y=%e\n", code, x, y);
173
174 x = .999;
175 code = OneDFunctionOptimize(&y, &x, 0.11, -4, 4, trialFn, 50, 1e-6, 0, 1);
176 fprintf(stderr, "code: %ld x=%e, y=%e\n", code, x, y);
177 exit(EXIT_SUCCESS);
178 }
179#endif
180
181 argc =
scanargs(&scArg, argc, argv);
182 if (argc < 3) {
183 fprintf(stderr, "%s%s", USAGE1, USAGE2);
184 exit(EXIT_FAILURE);
185 }
186 output = input = NULL;
187 flags = pipeFlags = excludes = truncate = pairFlags = 0;
188 indepQuantity = NULL;
189 depenQuantity = exclude = depenQuantityPair = NULL;
190 depenQuantities = 0;
191 noWarnings = 0;
192 fracRMSChangeLimit = 0.0;
193 fracFreqAccuracyLimit = 0.00001;
194 frequenciesDesired = 1;
195 maxFrequencies = 4;
196 freqCycleLimit = 100;
197 pairs = 0;
198
199 for (iArg = 1; iArg < argc; iArg++) {
200 if (scArg[iArg].arg_type == OPTION) {
201
202 switch (
match_string(scArg[iArg].list[0], option, N_OPTIONS, 0)) {
203 case SET_MAJOR_ORDER:
204 majorOrderFlag = 0;
205 scArg[iArg].n_items--;
206 if (scArg[iArg].n_items > 0 &&
207 (!
scanItemList(&majorOrderFlag, scArg[iArg].list + 1, &scArg[iArg].n_items, 0,
208 "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
209 "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL))) {
210 SDDS_Bomb(
"invalid -majorOrder syntax/values");
211 }
212 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
213 columnMajorOrder = 1;
214 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
215 columnMajorOrder = 0;
216 break;
217
218 case SET_TRUNCATE:
219 truncate = 1;
220 break;
221
222 case SET_PAIR:
223 if (depenQuantities)
224 SDDS_Bomb(
"Invalid -pair option, the depen-quantity is provided by -columns option already.");
225 if (scArg[iArg].n_items != 3)
227 depenQuantity =
SDDS_Realloc(depenQuantity,
sizeof(*depenQuantity) * (pairs + 1));
228 depenQuantityPair =
SDDS_Realloc(depenQuantityPair,
sizeof(*depenQuantityPair) * (pairs + 1));
229 depenQuantity[pairs] = scArg[iArg].list[1];
230 depenQuantityPair[pairs] = scArg[iArg].list[2];
231 pairs++;
232 break;
233
234 case SET_COLUMN:
235 if (indepQuantity)
236 SDDS_Bomb(
"only one -columns option may be given");
237 if (scArg[iArg].n_items < 2)
239 indepQuantity = scArg[iArg].list[1];
240 if (scArg[iArg].n_items >= 2) {
241 if (pairs)
242 SDDS_Bomb(
"Invalid -columns syntax, the depen-quantity is provided by -pair option already.");
243 depenQuantity =
tmalloc(
sizeof(*depenQuantity) * (depenQuantities = scArg[iArg].n_items - 2));
244 for (i = 0; i < depenQuantities; i++)
245 depenQuantity[i] = scArg[iArg].list[i + 2];
246 }
247 break;
248
249 case SET_PIPE:
250 if (!
processPipeOption(scArg[iArg].list + 1, scArg[iArg].n_items - 1, &pipeFlags))
252 break;
253
254 case SET_EXCLUDE:
255 if (scArg[iArg].n_items < 2)
257 moveToStringArray(&exclude, &excludes, scArg[iArg].list + 1, scArg[iArg].n_items - 1);
258 break;
259
260 case SET_NOWARNINGS:
261 noWarnings = 1;
262 break;
263
264 case SET_TERM_SEARCH:
265 items = scArg[iArg].n_items - 1;
266 flags &= ~(NAFF_RMS_CHANGE_LIMIT | NAFF_FREQS_DESIRED | NAFF_MAX_FREQUENCIES);
267 fracRMSChangeLimit = 0;
268 frequenciesDesired = 0;
269 maxFrequencies = 10;
270 if (!
scanItemList(&tmpFlags, scArg[iArg].list + 1, &items, 0,
271 "changelimit",
SDDS_DOUBLE, &fracRMSChangeLimit, 1, NAFF_RMS_CHANGE_LIMIT,
272 "maxfrequencies",
SDDS_LONG, &maxFrequencies, 1, NAFF_MAX_FREQUENCIES,
273 "frequencies",
SDDS_LONG, &frequenciesDesired, 1, NAFF_FREQS_DESIRED, NULL) ||
274 (tmpFlags & NAFF_RMS_CHANGE_LIMIT && tmpFlags & NAFF_FREQS_DESIRED) ||
275 maxFrequencies < 1) {
276 SDDS_Bomb(
"invalid -terminateSearch syntax");
277 }
278 flags |= tmpFlags;
279 if (frequenciesDesired)
280 maxFrequencies = frequenciesDesired;
281 break;
282
283 case SET_ITERATE_FREQ:
284 items = scArg[iArg].n_items - 1;
285 flags &= ~(NAFF_FREQ_CYCLE_LIMIT | NAFF_FREQ_ACCURACY_LIMIT);
286 if (!
scanItemList(&tmpFlags, scArg[iArg].list + 1, &items, 0,
287 "cyclelimit",
SDDS_LONG, &freqCycleLimit, 1, NAFF_FREQ_CYCLE_LIMIT,
288 "accuracylimit",
SDDS_DOUBLE, &fracFreqAccuracyLimit, 1, NAFF_FREQ_ACCURACY_LIMIT, NULL) ||
290 freqCycleLimit < 2) {
291 SDDS_Bomb(
"invalid -iterateFrequency syntax");
292 }
293 flags |= tmpFlags;
294 break;
295
296 default:
297 fprintf(stderr, "Error: unknown or ambiguous option: %s\n", scArg[iArg].list[0]);
298 exit(EXIT_FAILURE);
299 break;
300 }
301 } else {
302 if (!input)
303 input = scArg[iArg].list[0];
304 else if (!output)
305 output = scArg[iArg].list[0];
306 else
307 SDDS_Bomb(
"too many filenames provided");
308 }
309 }
310
312
313 if (!indepQuantity)
314 SDDS_Bomb(
"Supply the independent quantity name with the -columns option");
315
318
320 exit(EXIT_FAILURE);
321
322 excludes = appendToStringArray(&exclude, excludes, indepQuantity);
323 if (pairs) {
324 pairFlags = flags | NAFF_FREQ_FOUND;
325 depenQuantities = pairs;
326 }
327 if (!depenQuantities)
328 depenQuantities = appendToStringArray(&depenQuantity, depenQuantities, "*");
329 if (!pairs) {
330 if ((depenQuantities = expandColumnPairNames(&SDDSin, &depenQuantity, NULL, depenQuantities, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
332 SDDS_Bomb(
"No quantities selected to FFT");
333 }
334 }
335
336 if (!SetupNAFFOutput(&SDDSout, output, &SDDSin, indepQuantity, depenQuantities, depenQuantity,
337 &frequencyIndex, &litudeIndex, &phaseIndex, &significanceIndex,
338 depenQuantityPair, &litudeIndex1, &phaseIndex1, &significanceIndex1,
339 columnMajorOrder)) {
341 }
342
343 if (!(frequency =
SDDS_Malloc(
sizeof(*frequency) * maxFrequencies)) ||
344 !(amplitude =
SDDS_Malloc(
sizeof(*amplitude) * maxFrequencies)) ||
345 !(phase =
SDDS_Malloc(
sizeof(*phase) * maxFrequencies)) ||
346 !(significance =
SDDS_Malloc(
sizeof(*significance) * maxFrequencies))) {
348 }
349 if (pairs) {
350 if (!(amplitude1 =
SDDS_Malloc(
sizeof(*amplitude1) * maxFrequencies)) ||
351 !(phase1 =
SDDS_Malloc(
sizeof(*phase1) * maxFrequencies)) ||
352 !(significance1 =
SDDS_Malloc(
sizeof(*significance1) * maxFrequencies))) {
354 }
355 }
356
360 if (rows) {
361 int64_t primeRows;
362 rowsToUse = rows;
363 primeRows = greatestProductOfSmallPrimes(rows);
364 if (rows != primeRows) {
365 if (truncate)
366 rowsToUse = greatestProductOfSmallPrimes(rows);
368 fputs("Warning: Number of points has large prime factors.\n"
369 "This could take a very long time.\nConsider using the -truncate option.\n",
370 stderr);
371 }
375 }
378 for (i = 1; i < rowsToUse; i++)
379 if (tdata[i] <= tdata[i - 1])
380 SDDS_Bomb(
"Independent data is not monotonically increasing");
381 dt = (tdata[rowsToUse - 1] - tdata[0]) / (rowsToUse - 1.0);
382 t0 = tdata[0];
383 free(tdata);
384 tdata = NULL;
385 for (i = 0; i < depenQuantities; i++) {
388 for (j = 0; j < maxFrequencies; j++)
389 frequency[j] = amplitude[j] = phase[j] = significance[j] = -1;
390 PerformNAFF(frequency, amplitude, phase, significance, t0, dt, data, rowsToUse, flags,
391 fracRMSChangeLimit, maxFrequencies, freqCycleLimit, fracFreqAccuracyLimit, 0, 0);
392#ifdef DEBUG
393 fprintf(stderr, "Column %s: ", depenQuantity[i]);
394 fprintf(stderr, "f=%10.3e a=%10.3e p=%10.3e s=%10.3e\n", frequency[0], amplitude[0], phase[0], significance[0]);
395#endif
396 free(data);
397 data = NULL;
398 if (pairs) {
401 PerformNAFF(frequency, amplitude1, phase1, significance1, t0, dt, data, rowsToUse,
402 pairFlags, fracRMSChangeLimit, maxFrequencies, freqCycleLimit, fracFreqAccuracyLimit, 0, 0);
403
404 for (j = 0; j < maxFrequencies; j++)
405 if (frequency[j] != -1)
406 frequency[j] = adjustFrequencyHalfPlane(frequency[j], phase[j], phase1[j], dt);
407 free(data);
408 data = NULL;
409 }
410 if (!
SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, frequency, maxFrequencies, frequencyIndex[i]) ||
411 !
SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, amplitude, maxFrequencies, amplitudeIndex[i]) ||
412 !
SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, phase, maxFrequencies, phaseIndex[i]) ||
413 !
SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, significance, maxFrequencies, significanceIndex[i])) {
415 }
416 if (pairs) {
417 if (!
SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, amplitude1, maxFrequencies, amplitudeIndex1[i]) ||
418 !
SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, phase1, maxFrequencies, phaseIndex1[i]) ||
419 !
SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, significance1, maxFrequencies, significanceIndex1[i])) {
421 }
422 }
423 }
424 } else {
427 }
428 }
431 }
432
435 exit(EXIT_FAILURE);
436 }
439 exit(EXIT_FAILURE);
440 }
441 free(frequency);
442 free(amplitude);
443 free(phase);
444 free(significance);
445 if (pairs) {
446 free(amplitude1);
447 free(phase1);
448 free(significance1);
449 free(amplitudeIndex1);
450 free(phaseIndex1);
451 free(significanceIndex1);
452 free(depenQuantityPair);
453 }
454 free(depenQuantity);
455 free(frequencyIndex);
456 free(amplitudeIndex);
457 free(phaseIndex);
458 free(significanceIndex);
459 return EXIT_SUCCESS;
460}
int32_t SDDS_CopyParameters(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
int32_t SDDS_SetColumn(SDDS_DATASET *SDDS_dataset, int32_t mode, void *data, int64_t rows,...)
Sets the values for one data column in the current data table of an SDDS dataset.
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
int32_t SDDS_CheckColumn(SDDS_DATASET *SDDS_dataset, char *name, char *units, int32_t type, FILE *fp_message)
Checks if a column exists in the SDDS dataset with the specified name, units, and type.
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
void * SDDS_Malloc(size_t size)
Allocates memory of a specified size.
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
#define SDDS_ANY_NUMERIC_TYPE
Special identifier used by SDDS_Check*() routines to accept any numeric type.
#define SDDS_DOUBLE
Identifier for the double data type.
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
long bitsSet(unsigned long data)
Counts the number of set bits (1s) in the given data.
int64_t largest_prime_factor(int64_t number)
Find the largest prime factor of a number.
long match_string(char *string, char **option, long n_options, long mode)
Matches a given string against an array of option strings based on specified modes.
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
long processPipeOption(char **item, long items, unsigned long *flags)
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
long scanItemList(unsigned long *flags, char **item, long *items, unsigned long mode,...)
Scans a list of items and assigns values based on provided keywords and types.