86char *option[N_OPTIONS] = {
99 "Usage: sddsnaff [<inputfile>] [<outputfile>]\n"
100 " [-pipe=[input][,output]]\n"
101 " [-columns=<indep-variable>[,<depen-quantity>[,...]]]\n"
102 " [-pair=<column1>,<column2>]\n"
103 " [-exclude=<depen-quantity>[,...]]\n"
104 " [-terminateSearch={changeLimit=<fraction>[,maxFrequencies=<number>] | frequencies=<number>}]\n"
105 " [-iterateFrequency=[cycleLimit=<number>][,accuracyLimit=<fraction>]]\n"
108 " [-majorOrder=row|column]\n\n"
109 "Determines frequency components of signals using Laskar's NAFF method.\n"
110 "FFTs are involved in this process, hence some parameters refer to FFT configurations.\n\n"
112 " -pipe Use standard SDDS Toolkit pipe option for input and/or output.\n"
113 " -columns Specify the independent variable and dependent quantities to analyze.\n"
114 " <depen-quantity> entries may include wildcards.\n"
115 " -pair Specify a pair of columns for frequency and phase analysis.\n"
116 " Multiple -pair options can be provided.\n"
117 " -exclude Exclude specified quantities from analysis using wildcard patterns.\n"
118 " -terminateSearch Terminate the search based on RMS change limit or a specific number of frequencies.\n"
119 " -iterateFrequency Configure iteration parameters for frequency determination.\n"
120 " -truncate Truncate data to optimize FFT performance.\n"
121 " -noWarnings Suppress warning messages.\n"
122 " -majorOrder Specify output file's data order as row-major or column-major.\n\n";
125 "Program by Michael Borland. (" __DATE__
" " __TIME__
", SVN revision: " SVN_VERSION
")\n";
128 char *indepQuantity,
long depenQuantities,
char **depenQuantity,
long **frequencyIndex,
129 long **amplitudeIndex,
long **phaseIndex,
long **significanceIndex,
char **depenQuantityPair,
130 long **amplitudeIndex1,
long **phaseIndex1,
long **significanceIndex1,
short columnMajorOrder);
133double trialFn(
double x,
long *invalid) {
138 return (sin(x) / x + sqr(x - 0.5) * x * 0.5);
142int main(
int argc,
char **argv) {
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;
151 double *tdata, *data, t0, dt;
152 double fracRMSChangeLimit, fracFreqAccuracyLimit;
153 int32_t frequenciesDesired, maxFrequencies, freqCycleLimit;
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;
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);
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);
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);
181 argc =
scanargs(&scArg, argc, argv);
183 fprintf(stderr,
"%s%s", USAGE1, USAGE2);
186 output = input = NULL;
187 flags = pipeFlags = excludes = truncate = pairFlags = 0;
188 indepQuantity = NULL;
189 depenQuantity = exclude = depenQuantityPair = NULL;
192 fracRMSChangeLimit = 0.0;
193 fracFreqAccuracyLimit = 0.00001;
194 frequenciesDesired = 1;
196 freqCycleLimit = 100;
199 for (iArg = 1; iArg < argc; iArg++) {
200 if (scArg[iArg].arg_type == OPTION) {
202 switch (
match_string(scArg[iArg].list[0], option, N_OPTIONS, 0)) {
203 case SET_MAJOR_ORDER:
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");
212 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
213 columnMajorOrder = 1;
214 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
215 columnMajorOrder = 0;
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];
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) {
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];
250 if (!
processPipeOption(scArg[iArg].list + 1, scArg[iArg].n_items - 1, &pipeFlags))
255 if (scArg[iArg].n_items < 2)
257 moveToStringArray(&exclude, &excludes, scArg[iArg].list + 1, scArg[iArg].n_items - 1);
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;
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");
279 if (frequenciesDesired)
280 maxFrequencies = frequenciesDesired;
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");
297 fprintf(stderr,
"Error: unknown or ambiguous option: %s\n", scArg[iArg].list[0]);
303 input = scArg[iArg].list[0];
305 output = scArg[iArg].list[0];
307 SDDS_Bomb(
"too many filenames provided");
314 SDDS_Bomb(
"Supply the independent quantity name with the -columns option");
322 excludes = appendToStringArray(&exclude, excludes, indepQuantity);
324 pairFlags = flags | NAFF_FREQ_FOUND;
325 depenQuantities = pairs;
327 if (!depenQuantities)
328 depenQuantities = appendToStringArray(&depenQuantity, depenQuantities,
"*");
330 if ((depenQuantities = expandColumnPairNames(&SDDSin, &depenQuantity, NULL, depenQuantities, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
332 SDDS_Bomb(
"No quantities selected to FFT");
336 if (!SetupNAFFOutput(&SDDSout, output, &SDDSin, indepQuantity, depenQuantities, depenQuantity,
337 &frequencyIndex, &litudeIndex, &phaseIndex, &significanceIndex,
338 depenQuantityPair, &litudeIndex1, &phaseIndex1, &significanceIndex1,
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))) {
350 if (!(amplitude1 =
SDDS_Malloc(
sizeof(*amplitude1) * maxFrequencies)) ||
351 !(phase1 =
SDDS_Malloc(
sizeof(*phase1) * maxFrequencies)) ||
352 !(significance1 =
SDDS_Malloc(
sizeof(*significance1) * maxFrequencies))) {
363 primeRows = greatestProductOfSmallPrimes(rows);
364 if (rows != primeRows) {
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",
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);
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);
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]);
401 PerformNAFF(frequency, amplitude1, phase1, significance1, t0, dt, data, rowsToUse,
402 pairFlags, fracRMSChangeLimit, maxFrequencies, freqCycleLimit, fracFreqAccuracyLimit, 0, 0);
404 for (j = 0; j < maxFrequencies; j++)
405 if (frequency[j] != -1)
406 frequency[j] = adjustFrequencyHalfPlane(frequency[j], phase[j], phase1[j], dt);
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])) {
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])) {
449 free(amplitudeIndex1);
451 free(significanceIndex1);
452 free(depenQuantityPair);
455 free(frequencyIndex);
456 free(amplitudeIndex);
458 free(significanceIndex);
463 char *indepQuantity,
long depenQuantities,
char **depenQuantity,
long **frequencyIndex,
464 long **amplitudeIndex,
long **phaseIndex,
long **significanceIndex,
char **depenQuantityPair,
465 long **amplitudeIndex1,
long **phaseIndex1,
long **significanceIndex1,
short columnMajorOrder) {
466 char *freqUnits, *buffer, *ampUnits;
467 long i, maxBuffer, bufferNeeded;
469 if (!(*frequencyIndex =
SDDS_Malloc(
sizeof(**frequencyIndex) * depenQuantities)) ||
470 !(*amplitudeIndex =
SDDS_Malloc(
sizeof(**amplitudeIndex) * depenQuantities)) ||
471 !(*phaseIndex =
SDDS_Malloc(
sizeof(**phaseIndex) * depenQuantities)) ||
472 !(*significanceIndex =
SDDS_Malloc(
sizeof(**significanceIndex) * depenQuantities)) ||
473 !(buffer =
SDDS_Malloc(
sizeof(*buffer) * (maxBuffer = 1024)))) {
477 if (!(freqUnits = makeFrequencyUnits(SDDSin, indepQuantity)) ||
481 if (columnMajorOrder != -1)
482 SDDSout->layout.data_mode.column_major = columnMajorOrder;
484 SDDSout->layout.data_mode.column_major = SDDSin->layout.data_mode.column_major;
486 for (i = 0; i < depenQuantities; i++) {
487 if ((bufferNeeded = strlen(depenQuantity[i]) + 12) > maxBuffer &&
488 !(buffer =
SDDS_Realloc(buffer,
sizeof(*buffer) * bufferNeeded))) {
491 sprintf(buffer,
"%sFrequency", depenQuantity[i]);
496 sprintf(buffer,
"%sAmplitude", depenQuantity[i]);
499 sprintf(buffer,
"%sPhase", depenQuantity[i]);
502 sprintf(buffer,
"%sSignificance", depenQuantity[i]);
506 if (depenQuantityPair) {
507 if (!(*amplitudeIndex1 =
SDDS_Malloc(
sizeof(**amplitudeIndex1) * depenQuantities)) ||
508 !(*phaseIndex1 =
SDDS_Malloc(
sizeof(**phaseIndex1) * depenQuantities)) ||
509 !(*significanceIndex1 =
SDDS_Malloc(
sizeof(**significanceIndex1) * depenQuantities))) {
513 for (i = 0; i < depenQuantities; i++) {
514 if ((bufferNeeded = strlen(depenQuantityPair[i]) + 12) > maxBuffer &&
515 !(buffer =
SDDS_Realloc(buffer,
sizeof(*buffer) * bufferNeeded))) {
518 sprintf(buffer,
"%sAmplitude", depenQuantityPair[i]);
521 sprintf(buffer,
"%sPhase", depenQuantityPair[i]);
524 sprintf(buffer,
"%sSignificance", depenQuantityPair[i]);
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
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_GetColumnInformation(SDDS_DATASET *SDDS_dataset, char *field_name, void *memory, int32_t mode,...)
Retrieves information about a specified column in the SDDS dataset.
int32_t SDDS_InitializeOutput(SDDS_DATASET *SDDS_dataset, int32_t data_mode, int32_t lines_per_row, const char *description, const char *contents, const char *filename)
Initializes the SDDS output dataset.
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
int32_t SDDS_DefineColumn(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, int32_t field_length)
Defines a data column within the SDDS dataset.
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
int32_t SDDS_TransferAllParameterDefinitions(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source, uint32_t mode)
Transfers all parameter definitions from a source dataset to a target dataset.
void SDDS_SetError(char *error_text)
Records an error message in the SDDS error stack.
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_STRING
Identifier for the string data type.
#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.
Utility functions for SDDS dataset manipulation and string array operations.
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.