SDDS ToolKit Programs and Libraries for C and Python
Loading...
Searching...
No Matches
sddsfft.c
Go to the documentation of this file.
1/**
2 * @file sddsfft.c
3 * @brief Performs FFT (Fast Fourier Transform) on SDDS-formatted data files.
4 *
5 * @details
6 * This program processes SDDS data files, performing Fast Fourier Transforms (FFT) with extensive
7 * options for data manipulation, including normalization, windowing, PSD computation, and inverse
8 * transformations. It supports both complex and real input data and allows flexible output formats
9 * based on user-defined parameters.
10 *
11 * @section Usage
12 * ```
13 * sddsfft [<inputfile>] [<outputfile>]
14 * [-pipe=[input][,output]]
15 * [-columns=<indep-variable>[,<depen-quantity>[,...]]]
16 * [-complexInput[=unfolded|folded]]
17 * [-exclude=<depen-quantity>[,...]]
18 * [-window[={hanning|welch|parzen|hamming|flattop|gaussian|none}[,correct]]]
19 * [-sampleInterval=<number>]
20 * [-normalize]
21 * [-fullOutput[=unfolded|folded],unwrapLimit=<value>]
22 * [-psdOutput[=plain][,{integrated|rintegrated[=<cutoff>]}]]
23 * [-inverse]
24 * [-padwithzeroes[=exponent] | -truncate]
25 * [-suppressaverage]
26 * [-noWarnings]
27 * [-majorOrder=row|column]
28 * [-threads=<number>]
29 * ```
30 *
31 * @section Options
32 * | Required | Description |
33 * |--------------------|----------------------------------------------------------------------------|
34 * | `-columns` | Specify the independent variable and dependent quantities for FFT analysis. |
35 *
36 * | Optional | Description |
37 * |---------------------------------------|-----------------------------------------------------------------------------|
38 * | `-pipe` | Utilize the standard SDDS Toolkit pipe option for input and/or output. |
39 * | `-complexInput` | Specify complex input column handling. |
40 * | `-exclude` | List of wildcard patterns to exclude specific quantities from analysis. |
41 * | `-window` | Apply a windowing function before analysis. |
42 * | `-sampleInterval` | Define the interval for sampling input data points. |
43 * | `-normalize` | Normalize output to have a peak magnitude of 1. |
44 * | `-fullOutput` | Output the real and imaginary parts of the FFT. |
45 * | `-psdOutput` | Output Power Spectral Density (PSD) in various formats. |
46 * | `-inverse` | Perform an inverse Fourier transform. |
47 * | `-padwithzeroes` | Pad data with zeroes for FFT optimization. |
48 * | `-truncate` | Truncate data to nearest product of small primes for efficiency. |
49 * | `-suppressaverage` | Remove average value from data before FFT. |
50 * | `-noWarnings` | Suppress warning messages. |
51 * | `-majorOrder` | Specify row or column major order for output. |
52 * | `-threads` | Number of threads for per-column FFT preparation and result calculations. |
53 *
54 * @subsection Incompatibilities
55 * - `-truncate` is incompatible with:
56 * - `-padwithzeroes`
57 *
58 * @copyright
59 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
60 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
61 *
62 * @license
63 * This file is distributed under the terms of the Software License Agreement
64 * found in the file LICENSE included with this distribution.
65 *
66 * @authors
67 * - M. Borland
68 * - C. Saunders
69 * - R. Soliday
70 * - L. Emery
71 */
72
73#include "mdb.h"
74#include "SDDS.h"
75#include "scan.h"
76#include "fftpackC.h"
77#include "SDDSutils.h"
78#include <ctype.h>
79
80/* Enumeration for option types */
81enum option_type {
82 SET_WINDOW,
83 SET_NORMALIZE,
84 SET_PADWITHZEROES,
85 SET_TRUNCATE,
86 SET_SUPPRESSAVERAGE,
87 SET_SAMPLEINTERVAL,
88 SET_COLUMNS,
89 SET_FULLOUTPUT,
90 SET_PIPE,
91 SET_PSDOUTPUT,
92 SET_EXCLUDE,
93 SET_NOWARNINGS,
94 SET_COMPLEXINPUT,
95 SET_INVERSE,
96 SET_MAJOR_ORDER,
97 SET_THREADS,
98 N_OPTIONS
99};
100
101char *option[N_OPTIONS] = {
102 "window",
103 "normalize",
104 "padwithzeroes",
105 "truncate",
106 "suppressaverage",
107 "sampleinterval",
108 "columns",
109 "fulloutput",
110 "pipe",
111 "psdoutput",
112 "exclude",
113 "nowarnings",
114 "complexinput",
115 "inverse",
116 "majorOrder",
117 "threads",
118};
119
120#define FL_TRUNCATE 0x0001
121#define FL_PADWITHZEROES 0x0002
122#define FL_NORMALIZE 0x0004
123#define FL_SUPPRESSAVERAGE 0x0008
124#define FL_FULLOUTPUT 0x0010
125#define FL_MAKEFREQDATA 0x0020
126#define FL_PSDOUTPUT 0x0040
127#define FL_PSDINTEGOUTPUT 0x0080
128#define FL_PSDRINTEGOUTPUT 0x0100
129#define FL_FULLOUTPUT_FOLDED 0x0200
130#define FL_FULLOUTPUT_UNFOLDED 0x0400
131#define FL_COMPLEXINPUT_FOLDED 0x0800
132#define FL_COMPLEXINPUT_UNFOLDED 0x1000
133#define FL_UNWRAP_PHASE 0x2000
134
135#define WINDOW_HANNING 0
136#define WINDOW_WELCH 1
137#define WINDOW_PARZEN 2
138#define WINDOW_HAMMING 3
139#define WINDOW_FLATTOP 4
140#define WINDOW_GAUSSIAN 5
141#define WINDOW_NONE 6
142#define N_WINDOW_TYPES 7
143char *window_type[N_WINDOW_TYPES] = {
144 "hanning", "welch", "parzen", "hamming", "flattop", "gaussian", "none"};
145
146/* Improved and more readable usage message */
147static char *USAGE1 =
148 "Usage:\n"
149 " sddsfft [<inputfile>] [<outputfile>]\n"
150 " [-pipe=[input][,output]]\n"
151 " [-columns=<indep-variable>[,<depen-quantity>[,...]]]\n"
152 " [-complexInput[=unfolded|folded]]\n"
153 " [-exclude=<depen-quantity>[,...]]\n"
154 " [-window[={hanning|welch|parzen|hamming|flattop|gaussian|none}[,correct]]]\n"
155 " [-sampleInterval=<number>]\n"
156 " [-normalize]\n"
157 " [-fullOutput[=unfolded|folded],unwrapLimit=<value>]\n"
158 " [-psdOutput[=plain][,{integrated|rintegrated[=<cutoff>]}]]\n"
159 " [-inverse]\n"
160 " [-padwithzeroes[=exponent] | -truncate]\n"
161 " [-suppressaverage]\n"
162 " [-noWarnings]\n"
163 " [-majorOrder=row|column]\n"
164 " [-threads=<number>]\n\n";
165
166static char *USAGE2 =
167 "Options:\n"
168 " -pipe\n"
169 " Utilize the standard SDDS Toolkit pipe option for input and/or output.\n\n"
170 " -columns\n"
171 " Specify the independent variable and dependent quantities to Fourier analyze.\n"
172 " <depen-quantity> entries may include wildcards.\n\n"
173 " -complexInput\n"
174 " Indicate that input columns are complex, with names prefixed by Real and Imag.\n"
175 " Options:\n"
176 " folded (default): Input frequency space is folded.\n"
177 " unfolded : Input frequency space is unfolded and must include negative frequencies.\n"
178 " If omitted, the program checks the SpectrumFolded parameter in the input file.\n\n"
179 " -exclude\n"
180 " Provide a list of wildcard patterns to exclude specific quantities from analysis.\n\n"
181 " -window\n"
182 " Apply a windowing function to the data before analysis.\n"
183 " Available types:\n"
184 " hanning, welch, parzen, hamming, flattop, gaussian, none\n"
185 " Adding ',correct' applies a correction factor to preserve the integrated PSD.\n"
186 " Default: hanning.\n\n"
187 " -sampleInterval\n"
188 " Sample the input data points at the specified interval.\n\n"
189 " -normalize\n"
190 " Normalize the output to have a peak magnitude of 1.\n\n"
191 " -fullOutput\n"
192 " Output the real and imaginary parts of the FFT.\n"
193 " Options:\n"
194 " folded (default): Outputs the half FFT spectrum.\n"
195 " unfolded : Outputs the full FFT spectrum.\n"
196 " Adding ',unwrapLimit=<value>' computes and outputs the unwrapped phase where the relative magnitude exceeds the limit.\n\n"
197 " -psdOutput\n"
198 " Output the Power Spectral Density (PSD).\n"
199 " Options:\n"
200 " plain : Outputs the standard PSD.\n"
201 " integrated : Outputs the integrated PSD.\n"
202 " rintegrated=<cutoff> : Outputs the reverse-integrated PSD with an optional cutoff frequency.\n"
203 " Multiple options can be combined using commas.\n\n"
204 " -inverse\n"
205 " Perform an inverse Fourier transform. The output will always be an unfolded spectrum.\n"
206 " If combined with -fullOutput=folded, it overrides to -fullOutput=unfolded.\n\n"
207 " -padwithzeroes\n"
208 " Pad data with zeroes to optimize FFT performance.\n"
209 " Optionally specify an exponent to determine the padding size as 2^(original points * exponent).\n"
210 " -truncate\n"
211 " Truncate data to the nearest product of small prime numbers to reduce runtime.\n"
212 " Note: Only one of -padwithzeroes or -truncate can be used.\n\n"
213 " -suppressaverage\n"
214 " Remove the average value from the data before performing the FFT.\n\n"
215 " -noWarnings\n"
216 " Suppress all warning messages.\n\n"
217 " -majorOrder\n"
218 " Specify the output file's data order:\n"
219 " row : Row-major order.\n"
220 " column : Column-major order.\n\n"
221 " -threads\n"
222 " Number of threads for per-column FFT preparation and result calculations.\n\n"
223 "Program by Michael Borland. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n";
224
225/* Rest of the code remains unchanged */
226
227/* ... [The rest of your code continues here without changes] ... */
228
229int64_t greatestProductOfSmallPrimes(int64_t rows);
230long process_data(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, double *tdata, int64_t rows, int64_t rowsToUse,
231 char *depenQuantity, char *depenQuantity2, unsigned long flags, long windowType,
232 int64_t sampleInterval, long correctWindowEffects, long inverse, double rintegCutOffFreq, double unwrapLimit,
233 int threads);
234long create_fft_frequency_column(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *timeName, char *freqUnits, long inverse);
235
236long create_fft_columns(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *origName, char *indepName,
237 char *freqUnits, long full_output, unsigned long psd_output, long complexInput, long inverse, long unwrap_phase);
238long create_fft_parameters(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *indepName, char *freqUnits);
239char *makeFrequencyUnits(SDDS_DATASET *SDDSin, char *indepName);
240long expandComplexColumnPairNames(SDDS_DATASET *SDDSin, char **name, char ***realName, char ***imagName,
241 long names, char **excludeName, long excludeNames, long typeMode, long typeValue);
242
243int main(int argc, char **argv) {
244 int iArg;
245 char *freqUnits;
246 char *indepQuantity, **depenQuantity, **exclude, **realQuan = NULL, **imagQuan = NULL;
247 long depenQuantities, excludes;
248 char *input, *output;
249 long i, readCode, noWarnings, complexInput, inverse, spectrumFoldParExist = 0;
250 int64_t rows, rowsToUse, sampleInterval;
251 int32_t spectrumFolded = 0, page = 0;
252 unsigned long flags, pipeFlags, complexInputFlags = 0, fullOutputFlags = 0, majorOrderFlag;
253 long windowType = -1;
254 SCANNED_ARG *scanned;
255 SDDS_DATASET SDDSin, SDDSout;
256 double *tdata, rintegCutOffFreq, unwrapLimit = 0;
257 long padFactor, correctWindowEffects = 0;
258 short columnMajorOrder = -1;
259 int threads = 1;
260
261 SDDS_RegisterProgramName(argv[0]);
262 argc = scanargs(&scanned, argc, argv);
263 if (argc < 3 || argc > (3 + N_OPTIONS)) {
264 fprintf(stderr, "%s%s", USAGE1, USAGE2);
265 exit(EXIT_FAILURE);
266 /* bomb(NULL, USAGE); */
267 }
268 rintegCutOffFreq = 0;
269 output = input = NULL;
270 flags = pipeFlags = excludes = complexInput = inverse = 0;
271 sampleInterval = 1;
272 indepQuantity = NULL;
273 depenQuantity = exclude = NULL;
274 depenQuantities = 0;
275 noWarnings = 0;
276 padFactor = 0;
277 for (iArg = 1; iArg < argc; iArg++) {
278 if (scanned[iArg].arg_type == OPTION) {
279 /* process options here */
280 switch (match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
281 case SET_NORMALIZE:
282 flags |= FL_NORMALIZE;
283 break;
284 case SET_WINDOW:
285 if (scanned[iArg].n_items != 1) {
286 if ((i = match_string(scanned[iArg].list[1], window_type, N_WINDOW_TYPES, 0)) < 0)
287 SDDS_Bomb("unknown window type");
288 windowType = i;
289 if (scanned[iArg].n_items > 2) {
290 if (strncmp(scanned[iArg].list[2], "correct", strlen(scanned[iArg].list[2])) == 0)
291 correctWindowEffects = 1;
292 else
293 SDDS_Bomb("invalid -window syntax");
294 }
295 } else
296 windowType = 0;
297 break;
298 case SET_PADWITHZEROES:
299 flags |= FL_PADWITHZEROES;
300 if (scanned[iArg].n_items != 1) {
301 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%ld", &padFactor) != 1 || padFactor < 1)
302 SDDS_Bomb("invalid -padwithzeroes syntax");
303 }
304 break;
305 case SET_TRUNCATE:
306 flags |= FL_TRUNCATE;
307 break;
308 case SET_SUPPRESSAVERAGE:
309 flags |= FL_SUPPRESSAVERAGE;
310 break;
311 case SET_SAMPLEINTERVAL:
312 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%" SCNd64, &sampleInterval) != 1 || sampleInterval <= 0)
313 SDDS_Bomb("invalid -sampleinterval syntax");
314 break;
315 case SET_COLUMNS:
316 if (indepQuantity)
317 SDDS_Bomb("only one -columns option may be given");
318 if (scanned[iArg].n_items < 2)
319 SDDS_Bomb("invalid -columns syntax");
320 indepQuantity = scanned[iArg].list[1];
321 if (scanned[iArg].n_items >= 2) {
322 depenQuantity = tmalloc(sizeof(*depenQuantity) * (depenQuantities = scanned[iArg].n_items - 2));
323 for (i = 0; i < depenQuantities; i++)
324 depenQuantity[i] = scanned[iArg].list[i + 2];
325 }
326 break;
327 case SET_FULLOUTPUT:
328 flags |= FL_FULLOUTPUT;
329 if (scanned[iArg].n_items >= 2) {
330 scanned[iArg].n_items--;
331 if (!scanItemList(&fullOutputFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "folded", -1, NULL, 0, FL_FULLOUTPUT_FOLDED, "unfolded", -1, NULL, 0, FL_FULLOUTPUT_UNFOLDED, "unwrapLimit", SDDS_DOUBLE, &unwrapLimit, 0, FL_UNWRAP_PHASE, NULL))
332 SDDS_Bomb("Invalid -fullOutput syntax");
333 scanned[iArg].n_items++;
334 if (fullOutputFlags & FL_FULLOUTPUT_UNFOLDED)
335 flags |= FL_FULLOUTPUT_UNFOLDED;
336 else
337 flags |= FL_FULLOUTPUT_FOLDED;
338 if (fullOutputFlags & FL_UNWRAP_PHASE)
339 flags |= FL_UNWRAP_PHASE;
340 }
341 break;
342 case SET_PIPE:
343 if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
344 SDDS_Bomb("invalid -pipe syntax");
345 break;
346 case SET_PSDOUTPUT:
347 if (scanned[iArg].n_items -= 1) {
348 unsigned long tmpFlags;
349 if (strchr(scanned[iArg].list[1], '=') <= 0) {
350 if (!scanItemList(&tmpFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "integrated", -1, NULL, 0, FL_PSDINTEGOUTPUT, "rintegrated", -1, NULL, 0, FL_PSDRINTEGOUTPUT, "plain", -1, NULL, 0, FL_PSDOUTPUT, NULL))
351 SDDS_Bomb("invalid -psdOutput syntax");
352 } else {
353 if (!scanItemList(&tmpFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "integrated", -1, NULL, 0, FL_PSDINTEGOUTPUT, "rintegrated", SDDS_DOUBLE, &rintegCutOffFreq, 0, FL_PSDRINTEGOUTPUT, "plain", -1, NULL, 0, FL_PSDOUTPUT, NULL))
354 SDDS_Bomb("invalid -psdOutput syntax");
355 }
356 flags |= tmpFlags;
357 } else
358 flags |= FL_PSDOUTPUT;
359 if (flags & FL_PSDINTEGOUTPUT && flags & FL_PSDRINTEGOUTPUT)
360 SDDS_Bomb("invalid -psdOutput syntax: give only one of integrated or rintegrated");
361 break;
362 case SET_EXCLUDE:
363 if (scanned[iArg].n_items < 2)
364 SDDS_Bomb("invalid -exclude syntax");
365 moveToStringArray(&exclude, &excludes, scanned[iArg].list + 1, scanned[iArg].n_items - 1);
366 break;
367 case SET_NOWARNINGS:
368 noWarnings = 1;
369 break;
370 case SET_COMPLEXINPUT:
371 complexInput = 1;
372 if (scanned[iArg].n_items == 2) {
373 scanned[iArg].n_items--;
374 if (!scanItemList(&complexInputFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "folded", -1, NULL, 0, FL_COMPLEXINPUT_FOLDED, "unfolded", -1, NULL, 0, FL_COMPLEXINPUT_UNFOLDED, NULL))
375 SDDS_Bomb("Invalid -complexInput syntax");
376 scanned[iArg].n_items++;
377 }
378 break;
379 case SET_INVERSE:
380 inverse = 1;
381 break;
382 case SET_MAJOR_ORDER:
383 majorOrderFlag = 0;
384 scanned[iArg].n_items--;
385 if (scanned[iArg].n_items > 0 && (!scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER, "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
386 SDDS_Bomb("invalid -majorOrder syntax/values");
387 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
388 columnMajorOrder = 1;
389 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
390 columnMajorOrder = 0;
391 break;
392 case SET_THREADS:
393 if (scanned[iArg].n_items != 2 ||
394 sscanf(scanned[iArg].list[1], "%d", &threads) != 1 || threads < 1)
395 SDDS_Bomb("invalid -threads syntax");
396 break;
397 default:
398 fprintf(stderr, "error: unknown/ambiguous option: %s\n", scanned[iArg].list[0]);
399 exit(EXIT_FAILURE);
400 break;
401 }
402 } else {
403 if (!input)
404 input = scanned[iArg].list[0];
405 else if (!output)
406 output = scanned[iArg].list[0];
407 else
408 SDDS_Bomb("too many filenames seen");
409 }
410 }
411 if (!complexInput) {
412 if (!noWarnings && inverse)
413 fprintf(stderr, "Warning: the inverse option is ignored since it only works with -complexInput.\n");
414 inverse = 0;
415 }
416 if (!noWarnings && inverse && flags & FL_FULLOUTPUT_FOLDED)
417 fprintf(stderr, "Warning: the -inverse -fullOutput=folded will be changed to -inverse -fullOutput=unfolded.\n");
418
419 processFilenames("sddsfft", &input, &output, pipeFlags, 0, NULL);
420
421 if (!indepQuantity)
422 SDDS_Bomb("Supply the independent quantity name with the -columns option.");
423
424 if (flags & FL_TRUNCATE && flags & FL_PADWITHZEROES)
425 SDDS_Bomb("Specify only one of -padwithzeroes and -truncate.");
426
427 if (!SDDS_InitializeInput(&SDDSin, input))
428 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
429
430 if (SDDS_CheckColumn(&SDDSin, indepQuantity, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OKAY)
431 exit(EXIT_FAILURE);
432
433 excludes = appendToStringArray(&exclude, excludes, indepQuantity);
434 if (!depenQuantities)
435 depenQuantities = appendToStringArray(&depenQuantity, depenQuantities, "*");
436
437 if (!complexInput) {
438 if ((depenQuantities = expandColumnPairNames(&SDDSin, &depenQuantity, NULL, depenQuantities, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
439 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
440 SDDS_Bomb("No quantities selected to FFT.");
441 }
442 } else {
443 if ((depenQuantities = expandComplexColumnPairNames(&SDDSin, depenQuantity, &realQuan, &imagQuan, depenQuantities, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
444 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
445 SDDS_Bomb("No quantities selected to FFT.");
446 }
447 }
448
449#if 0
450 fprintf(stderr, "%ld dependent quantities:\n", depenQuantities);
451 for (i = 0; i < depenQuantities; i++)
452 fprintf(stderr, " %s\n", depenQuantity[i]);
453#endif
454
455 if (!(freqUnits = makeFrequencyUnits(&SDDSin, indepQuantity)) ||
456 !SDDS_InitializeOutput(&SDDSout, SDDS_BINARY, 0, NULL, "sddsfft output", output) ||
457 !create_fft_frequency_column(&SDDSout, &SDDSin, indepQuantity, freqUnits, inverse) ||
458 SDDS_DefineParameter(&SDDSout, "fftFrequencies", NULL, NULL, NULL, NULL, SDDS_LONG, NULL) < 0 ||
459 SDDS_DefineParameter(&SDDSout, "fftFrequencySpacing", "$gD$rf", freqUnits, NULL, NULL, SDDS_DOUBLE, NULL) < 0)
460 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
461 if (columnMajorOrder != -1)
462 SDDSout.layout.data_mode.column_major = columnMajorOrder;
463 else
464 SDDSout.layout.data_mode.column_major = SDDSin.layout.data_mode.column_major;
465
466 if (flags & FL_FULLOUTPUT && SDDS_DefineParameter(&SDDSout, "SpectrumFolded", NULL, NULL, NULL, NULL, SDDS_LONG, NULL) < 0)
467 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
468 if (complexInput) {
469 if (!complexInputFlags) {
470 if (SDDS_CheckParameter(&SDDSin, "SpectrumFolded", NULL, SDDS_LONG, NULL) == SDDS_CHECK_OK)
471 spectrumFoldParExist = 1;
472 } else if (complexInputFlags & FL_COMPLEXINPUT_UNFOLDED)
473 flags |= FL_COMPLEXINPUT_UNFOLDED;
474 else
475 flags |= FL_COMPLEXINPUT_FOLDED;
476 }
477 for (i = 0; i < depenQuantities; i++) {
478 if (!complexInput)
479 create_fft_columns(&SDDSout, &SDDSin, depenQuantity[i], indepQuantity, freqUnits,
480 flags & FL_FULLOUTPUT,
481 flags & (FL_PSDOUTPUT + FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT),
482 0, inverse, flags & FL_UNWRAP_PHASE);
483 else
484 create_fft_columns(&SDDSout, &SDDSin, realQuan[i], indepQuantity, freqUnits,
485 flags & FL_FULLOUTPUT,
486 flags & (FL_PSDOUTPUT + FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT),
487 1, inverse, flags & FL_UNWRAP_PHASE);
488 }
489
490 if (!SDDS_TransferAllParameterDefinitions(&SDDSout, &SDDSin, SDDS_TRANSFER_KEEPOLD) ||
491 !SDDS_WriteLayout(&SDDSout))
492 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
493
494 while ((readCode = SDDS_ReadPage(&SDDSin)) > 0) {
495 page++;
496 if ((rows = SDDS_CountRowsOfInterest(&SDDSin)) < 0)
497 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
498 if (page == 1 && spectrumFoldParExist) {
499 if (!SDDS_GetParameterAsLong(&SDDSin, "SpectrumFolded", &spectrumFolded))
500 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
501 if (spectrumFolded)
502 flags |= FL_COMPLEXINPUT_FOLDED;
503 else
504 flags |= FL_COMPLEXINPUT_UNFOLDED;
505 }
506 if (rows) {
507 int64_t primeRows, pow2Rows;
508 if (rows < 2)
509 SDDS_Bomb("sddsfft requires at least two rows on each nonempty page");
510 rowsToUse = rows;
511 primeRows = greatestProductOfSmallPrimes(rows);
512 if (rows != primeRows || padFactor) {
513 if (flags & FL_PADWITHZEROES) {
514 pow2Rows = ipow(2., ((int64_t)(log((double)rows) / log(2.0F))) + (padFactor ? padFactor : 1.0));
515 if ((primeRows = greatestProductOfSmallPrimes(pow2Rows)) > rows)
516 rowsToUse = primeRows;
517 else
518 rowsToUse = pow2Rows;
519 } else if (flags & FL_TRUNCATE)
520 rowsToUse = greatestProductOfSmallPrimes(rows);
521 else if (largest_prime_factor(rows) > 1000 && !noWarnings)
522 fputs("Warning: number of points has large prime factors.\nThis could take a very long time.\nConsider using the -truncate option.\n", stderr);
523 }
524 if (!SDDS_StartPage(&SDDSout, 2 * rowsToUse + 2) || !SDDS_CopyParameters(&SDDSout, &SDDSin))
525 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
526 if (!(tdata = SDDS_GetColumnInDoubles(&SDDSin, indepQuantity)))
527 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
528 for (i = 0; i < depenQuantities; i++)
529 if (!process_data(&SDDSout, &SDDSin, tdata, rows, rowsToUse,
530 complexInput ? realQuan[i] : depenQuantity[i],
531 complexInput ? imagQuan[i] : NULL,
532 flags | (i == 0 ? FL_MAKEFREQDATA : 0),
533 windowType, sampleInterval, correctWindowEffects, inverse,
534 rintegCutOffFreq, unwrapLimit, threads)) {
535 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
536 exit(EXIT_FAILURE);
537 }
538 free(tdata);
539 } else {
540 if (!SDDS_StartPage(&SDDSout, 0) || !SDDS_CopyParameters(&SDDSout, &SDDSin))
541 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
542 }
543 if (!SDDS_WritePage(&SDDSout))
544 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
545 }
546
547 if (!SDDS_Terminate(&SDDSin)) {
548 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
549 exit(EXIT_FAILURE);
550 }
551 if (!SDDS_Terminate(&SDDSout)) {
552 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
553 exit(EXIT_FAILURE);
554 }
555
556 return EXIT_SUCCESS;
557}
558
559static long psdOffset, argOffset, realOffset, imagOffset, fftOffset = -1, psdIntOffset, psdIntPowerOffset, unwrappedArgOffset = -1;
560
561long process_data(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, double *tdata, int64_t rows,
562 int64_t rowsToUse, char *depenQuantity, char *imagQuantity, unsigned long flags, long windowType,
563 int64_t sampleInterval, long correctWindowEffects, long inverse, double rintegCutOffFreq, double unwrapLimit,
564 int threads) {
565 long offset, index, unfold = 0;
566 int64_t n_freq, i, fftrows = 0;
567 double r, r1, r2, length, factor, df, min, max, delta;
568 double *real, *imag, *magData, *arg = NULL, *real_imag, *data, *psd = NULL, *psdInteg = NULL, *psdIntegPower = NULL, *unwrapArg = NULL, phase_correction = 0;
569 double t0;
570 char s[256];
571 double *fdata, *imagData = NULL;
572 double *tDataStore = NULL;
573 double windowCorrectionFactor = 0;
574
575 if (!(data = SDDS_GetColumnInDoubles(SDDSin, depenQuantity)))
576 return 0;
577 if (imagQuantity && !(imagData = SDDS_GetColumnInDoubles(SDDSin, imagQuantity)))
578 return 0;
579 if (flags & FL_SUPPRESSAVERAGE) {
580 compute_average(&r, data, rows);
581#pragma omp parallel for if (threads > 1) num_threads(threads)
582 for (i = 0; i < rows; i++)
583 data[i] -= r;
584 if (imagData) {
585 compute_average(&r, imagData, rows);
586#pragma omp parallel for if (threads > 1) num_threads(threads)
587 for (i = 0; i < rows; i++)
588 imagData[i] -= r;
589 }
590 }
591 if (rows < rowsToUse) {
592 /* pad with zeroes */
593 tDataStore = tmalloc(sizeof(*tDataStore) * rowsToUse);
594 memcpy((char *)tDataStore, (char *)tdata, rows * sizeof(*tdata));
595 if (!(data = SDDS_Realloc(data, sizeof(*data) * rowsToUse)))
596 SDDS_Bomb("memory allocation failure");
597 if (imagData && !(imagData = SDDS_Realloc(imagData, sizeof(*imagData) * rowsToUse)))
598 SDDS_Bomb("memory allocation failure");
599 length = tdata[rows - 1] - tdata[0];
600 for (i = rows; i < rowsToUse; i++) {
601 tDataStore[i] = tDataStore[i - 1] + length / ((double)rows - 1);
602 data[i] = 0;
603 }
604 if (imagData)
605 for (i = rows; i < rowsToUse; i++)
606 imagData[i] = 0;
607 tdata = tDataStore;
608 }
609 rows = rowsToUse; /* results in truncation if rows>rowsToUse */
610 windowCorrectionFactor = 0;
611 switch (windowType) {
612 case WINDOW_HANNING:
613 r = PIx2 / (rows - 1);
614#pragma omp parallel for private(factor) reduction(+ : windowCorrectionFactor) if (threads > 1) num_threads(threads)
615 for (i = 0; i < rows; i++) {
616 factor = (1 - cos(i * r)) / 2;
617 data[i] *= factor;
618 windowCorrectionFactor += sqr(factor);
619 if (imagData)
620 imagData[i] *= factor;
621 }
622 break;
623 case WINDOW_HAMMING:
624 r = PIx2 / (rows - 1);
625#pragma omp parallel for private(factor) reduction(+ : windowCorrectionFactor) if (threads > 1) num_threads(threads)
626 for (i = 0; i < rows; i++) {
627 factor = 0.54 - 0.46 * cos(i * r);
628 data[i] *= factor;
629 windowCorrectionFactor += sqr(factor);
630 if (imagData)
631 imagData[i] *= factor;
632 }
633 if (imagData)
634#pragma omp parallel for if (threads > 1) num_threads(threads)
635 for (i = 0; i < rows; i++)
636 imagData[i] *= (1 - cos(i * r)) / 2;
637 break;
638 case WINDOW_WELCH:
639 r1 = (rows - 1) / 2.0;
640 r2 = sqr((rows + 1) / 2.0);
641#pragma omp parallel for private(factor) reduction(+ : windowCorrectionFactor) if (threads > 1) num_threads(threads)
642 for (i = 0; i < rows; i++) {
643 factor = 1 - sqr(i - r1) / r2;
644 data[i] *= factor;
645 windowCorrectionFactor += sqr(factor);
646 if (imagData)
647 imagData[i] *= factor;
648 }
649 break;
650 case WINDOW_PARZEN:
651 r = (rows - 1) / 2.0;
652#pragma omp parallel for private(factor) reduction(+ : windowCorrectionFactor) if (threads > 1) num_threads(threads)
653 for (i = 0; i < rows; i++) {
654 factor = 1 - FABS((i - r) / r);
655 data[i] *= factor;
656 windowCorrectionFactor += sqr(factor);
657 if (imagData)
658 imagData[i] *= factor;
659 }
660 break;
661 case WINDOW_FLATTOP:
662#pragma omp parallel for private(r, factor) reduction(+ : windowCorrectionFactor) if (threads > 1) num_threads(threads)
663 for (i = 0; i < rows; i++) {
664 r = i * PIx2 / (rows - 1);
665 factor = 1 - 1.93 * cos(r) + 1.29 * cos(2 * r) - 0.388 * cos(3 * r) + 0.032 * cos(4 * r);
666 data[i] *= factor;
667 windowCorrectionFactor += sqr(factor);
668 if (imagData)
669 imagData[i] *= factor;
670 }
671 break;
672 case WINDOW_GAUSSIAN:
673#pragma omp parallel for private(r, factor) reduction(+ : windowCorrectionFactor) if (threads > 1) num_threads(threads)
674 for (i = 0; i < rows; i++) {
675 r = sqr((i - (rows - 1) / 2.) / (0.4 * (rows - 1) / 2.)) / 2;
676 factor = exp(-r);
677 data[i] *= factor;
678 windowCorrectionFactor += sqr(factor);
679 if (imagData)
680 imagData[i] *= factor;
681 }
682 break;
683 case WINDOW_NONE:
684 default:
685 windowCorrectionFactor = 1;
686 break;
687 }
688
689 if (correctWindowEffects) {
690 /* Add correction factor to make the integrated PSD come out right. */
691 windowCorrectionFactor = 1 / sqrt(windowCorrectionFactor / rows);
692#pragma omp parallel for if (threads > 1) num_threads(threads)
693 for (i = 0; i < rows; i++)
694 data[i] *= windowCorrectionFactor;
695 if (imagData) {
696#pragma omp parallel for if (threads > 1) num_threads(threads)
697 for (i = 0; i < rows; i++)
698 imagData[i] *= windowCorrectionFactor;
699 }
700 }
701 if (imagData && flags & FL_COMPLEXINPUT_FOLDED) {
702 double min, max, max1;
703 data = SDDS_Realloc(data, sizeof(*data) * rows * 2);
704 imagData = SDDS_Realloc(imagData, sizeof(*data) * rows * 2);
705 find_min_max(&min, &max, data, rows);
706 if (fabs(min) > fabs(max))
707 max1 = fabs(min);
708 else
709 max1 = fabs(max);
710 find_min_max(&min, &max, imagData, rows);
711 if (fabs(min) > max1)
712 max1 = fabs(min);
713 if (fabs(max) > max1)
714 max1 = fabs(max);
715 if (fabs(imagData[rows - 1]) / max1 < 1.0e-15) {
716 fftrows = 2 * (rows - 1);
717 for (i = 1; i < rows - 1; i++) {
718 data[i] = data[i] / 2.0;
719 imagData[i] = imagData[i] / 2.0;
720 }
721 for (i = 1; i < rows - 1; i++) {
722 data[rows - 1 + i] = data[rows - 1 - i];
723 imagData[rows - 1 + i] = -imagData[rows - 1 - i];
724 }
725 length = (tdata[rows - 1] - tdata[0]) * 2.0;
726 } else {
727 fftrows = 2 * (rows - 1) + 1;
728 for (i = 1; i < rows; i++) {
729 data[i] = data[i] / 2.0;
730 imagData[i] = imagData[i] / 2.0;
731 }
732 for (i = 0; i < rows - 1; i++) {
733 data[rows + i] = data[rows - 1 - i];
734 imagData[rows + i] = -imagData[rows - 1 - i];
735 }
736 length = ((double)fftrows) * (tdata[rows - 1] - tdata[0]) / ((double)fftrows - 1.0) * 2;
737 }
738 } else {
739 fftrows = rows;
740 length = ((double)rows) * (tdata[rows - 1] - tdata[0]) / ((double)rows - 1.0);
741 }
742 /* compute FFT */
743
744 real_imag = tmalloc(sizeof(double) * (2 * fftrows + 2));
745#pragma omp parallel for if (threads > 1) num_threads(threads)
746 for (i = 0; i < fftrows; i++) {
747 real_imag[2 * i] = data[i];
748 if (imagData)
749 real_imag[2 * i + 1] = imagData[i];
750 else
751 real_imag[2 * i + 1] = 0;
752 }
753 if (!inverse) {
754 complexFFT(real_imag, fftrows, 0);
755 if (flags & FL_FULLOUTPUT_UNFOLDED) {
756 n_freq = fftrows;
757 unfold = 1;
758 } else if (flags & FL_FULLOUTPUT_FOLDED)
759 n_freq = fftrows / 2 + 1;
760 else if (!imagData)
761 n_freq = fftrows / 2 + 1;
762 else
763 n_freq = fftrows + 1;
764 } else {
765 complexFFT(real_imag, fftrows, INVERSE_FFT);
766 n_freq = fftrows;
767 }
768 /* calculate factor for converting k to f or omega */
769 /* length is assumed the length of period, not the total length of the data file */
770
771 /* length = ((double)rows)*(tdata[rows-1]-tdata[0])/((double)rows-1.0); */
772 t0 = tdata[0];
773 df = factor = 1.0 / length;
774
775 /* convert into amplitudes and frequencies, adding phase factor for t[0]!=0 */
776 real = tmalloc(sizeof(double) * n_freq);
777 imag = tmalloc(sizeof(double) * n_freq);
778 fdata = tmalloc(sizeof(double) * n_freq);
779 magData = tmalloc(sizeof(double) * n_freq);
780 if (flags & FL_PSDOUTPUT || flags & (FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT)) {
781 psd = tmalloc(sizeof(*psd) * n_freq);
782 if (flags & (FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT)) {
783 psdInteg = tmalloc(sizeof(*psdInteg) * n_freq);
784 psdIntegPower = tmalloc(sizeof(*psdIntegPower) * n_freq);
785 }
786 }
787
788#pragma omp parallel for if (threads > 1) num_threads(threads)
789 for (i = 0; i < n_freq; i++) {
790 double dtf_real, dtf_imag;
791 fdata[i] = i * df;
792 dtf_real = cos(-2 * PI * fdata[i] * t0);
793 dtf_imag = sin(-2 * PI * fdata[i] * t0);
794 if (psd)
795 psd[i] = (sqr(real_imag[2 * i]) + sqr(real_imag[2 * i + 1])) / df;
796 if (!imagData && i != 0 && !(i == (n_freq - 1) && rows % 2 == 0)) {
797 /* This is not the DC or Nyquist term, so
798 multiply by 2 to account for amplitude in
799 negative frequencies
800 */
801 if (!unfold) {
802 real_imag[2 * i] *= 2;
803 real_imag[2 * i + 1] *= 2;
804 }
805 if (psd)
806 psd[i] *= 2; /* 2 really is what I want--not 4 */
807 }
808 real[i] = real_imag[2 * i] * dtf_real - real_imag[2 * i + 1] * dtf_imag;
809 imag[i] = real_imag[2 * i + 1] * dtf_real + real_imag[2 * i] * dtf_imag;
810 magData[i] = sqrt(sqr(real[i]) + sqr(imag[i]));
811 }
812
813 if (psdInteg) {
814 if (flags & FL_PSDINTEGOUTPUT) {
815 psdIntegPower[0] = 0;
816 for (i = 1; i < n_freq; i++)
817 psdIntegPower[i] = psdIntegPower[i - 1] + (psd[i - 1] + psd[i]) * df / 2.;
818 for (i = 0; i < n_freq; i++)
819 psdInteg[i] = sqrt(psdIntegPower[i]);
820 } else {
821 psdIntegPower[n_freq - 1] = 0;
822 for (i = n_freq - 2; i >= 0; i--) {
823 if (rintegCutOffFreq == 0 || fdata[i] <= rintegCutOffFreq)
824 psdIntegPower[i] = psdIntegPower[i + 1] + (psd[i + 1] + psd[i]) * df / 2.;
825 }
826 for (i = 0; i < n_freq; i++)
827 psdInteg[i] = sqrt(psdIntegPower[i]);
828 }
829 }
830
831 if (flags & FL_FULLOUTPUT) {
832 arg = tmalloc(sizeof(*arg) * n_freq);
833#pragma omp parallel for if (threads > 1) num_threads(threads)
834 for (i = 0; i < n_freq; i++) {
835 if (real[i] || imag[i])
836 arg[i] = 180.0 / PI * atan2(imag[i], real[i]);
837 else
838 arg[i] = 0;
839 }
840 }
841 if (flags & FL_UNWRAP_PHASE) {
842 find_min_max(&min, &max, magData, n_freq);
843 unwrapArg = tmalloc(sizeof(*unwrapArg) * n_freq);
844 phase_correction = 0;
845 for (i = 0; i < n_freq; i++) {
846 if (i && magData[i] / max > unwrapLimit) {
847 delta = arg[i] - arg[i - 1];
848 if (delta < -180.0)
849 phase_correction += 360.0;
850 else if (delta > 180.0)
851 phase_correction -= 360.0;
852 }
853 unwrapArg[i] = arg[i] + phase_correction;
854 }
855 }
856
857 if (flags & FL_NORMALIZE) {
858 factor = -DBL_MAX;
859 for (i = 0; i < n_freq; i++)
860 if (magData[i] > factor)
861 factor = magData[i];
862 if (factor != -DBL_MAX)
863 for (i = 0; i < n_freq; i++) {
864 real[i] /= factor;
865 imag[i] /= factor;
866 magData[i] /= factor;
867 }
868 }
869 if (!inverse)
870 sprintf(s, "FFT%s", depenQuantity + (imagData ? 4 : 0));
871 else {
872 if (strncmp(depenQuantity, "FFT", 3) == 0)
873 sprintf(s, "%s", depenQuantity + 3);
874 else if (strncmp(depenQuantity, "RealFFT", 7) == 0)
875 sprintf(s, "%s", depenQuantity + 7);
876 else
877 sprintf(s, "%s", depenQuantity);
878 }
879
880 if ((index = SDDS_GetColumnIndex(SDDSout, s)) < 0)
881 return 0;
882
883 if (flags & FL_SUPPRESSAVERAGE) {
884 n_freq -= 1;
885 offset = 1;
886 } else
887 offset = 0;
888
889 if ((flags & FL_MAKEFREQDATA &&
890 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, fdata + offset, n_freq, 0)) ||
891 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, magData + offset, n_freq, index + fftOffset) ||
892 (flags & FL_FULLOUTPUT &&
893 (!SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, real + offset, n_freq, index + realOffset) ||
894 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, imag + offset, n_freq, index + imagOffset) ||
895 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, arg + offset, n_freq, index + argOffset))) ||
896 (flags & FL_PSDOUTPUT &&
897 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, psd + offset, n_freq, index + psdOffset)) ||
898 (flags & (FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT) && (!SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, psdInteg + offset, n_freq, index + psdIntOffset) ||
899 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, psdIntegPower + offset, n_freq, index + psdIntPowerOffset))) ||
900 (flags & FL_UNWRAP_PHASE && !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, unwrapArg + offset, n_freq, index + unwrappedArgOffset)))
901 return 0;
902 if (sampleInterval > 0) {
903 int32_t *sample_row_flag;
904 sample_row_flag = calloc(sizeof(*sample_row_flag), n_freq);
905 for (i = 0; i < n_freq; i += sampleInterval)
906 sample_row_flag[i] = 1;
907 if (!SDDS_AssertRowFlags(SDDSout, SDDS_FLAG_ARRAY, sample_row_flag, n_freq)) {
908 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
909 return 0;
910 }
911 free(sample_row_flag);
912 }
913 if (!SDDS_SetParameters(SDDSout, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, "fftFrequencies", n_freq, "fftFrequencySpacing", df, NULL))
914 return 0;
915 if (flags & FL_FULLOUTPUT && !SDDS_SetParameters(SDDSout, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, "SpectrumFolded", flags & FL_FULLOUTPUT_UNFOLDED ? 0 : 1, NULL))
916 return 0;
917 free(data);
918 free(magData);
919 if (imagData)
920 free(imagData);
921 if (tDataStore)
922 free(tDataStore);
923 if (arg)
924 free(arg);
925 free(real_imag);
926 free(real);
927 free(imag);
928 free(fdata);
929 if (psd)
930 free(psd);
931 if (psdInteg)
932 free(psdInteg);
933 if (psdIntegPower)
934 free(psdIntegPower);
935 if (unwrapArg)
936 free(unwrapArg);
937 return 1;
938}
939
940long create_fft_frequency_column(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *timeName, char *freqUnits, long inverse) {
941 char s[SDDS_MAXLINE];
942 char *timeSymbol;
943 char *description;
944
945 if (SDDS_GetColumnInformation(SDDSin, "symbol", &timeSymbol, SDDS_GET_BY_NAME, timeName) != SDDS_STRING)
946 return 0;
947 if (!timeSymbol || SDDS_StringIsBlank(timeSymbol))
948 SDDS_CopyString(&timeSymbol, timeName);
949
950 sprintf(s, "Frequency for %s", timeSymbol);
951 SDDS_CopyString(&description, s);
952 if (!inverse) {
953 if (SDDS_DefineColumn(SDDSout, "f", NULL, freqUnits, description, NULL, SDDS_DOUBLE, 0) < 0) {
954 free(timeSymbol);
955 free(description);
956 return 0;
957 }
958 } else {
959 sprintf(s, "inverse for %s", timeSymbol);
960 SDDS_CopyString(&description, s);
961 if (SDDS_DefineColumn(SDDSout, "t", NULL, freqUnits, description, NULL, SDDS_DOUBLE, 0) < 0) {
962 free(timeSymbol);
963 free(description);
964 return 0;
965 }
966 }
967 free(timeSymbol);
968 free(description);
969 return 1;
970}
971
972long create_fft_columns(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *origName, char *indepName, char *freqUnits, long full_output, unsigned long psd_output, long complexInput, long inverse, long unwrap_phase) {
973 char s[SDDS_MAXLINE];
974 char *origUnits, *origSymbol;
975 char *description, *name, *symbol, *units;
976 long index0, index1;
977 long offset = 0;
978
979 if (complexInput)
980 offset = 4;
981 if (SDDS_GetColumnInformation(SDDSin, "units", &origUnits, SDDS_GET_BY_NAME, origName) != SDDS_STRING ||
982 SDDS_GetColumnInformation(SDDSin, "symbol", &origSymbol, SDDS_GET_BY_NAME, origName) != SDDS_STRING)
983 return 0;
984 if (!inverse)
985 sprintf(s, "FFT%s", origName + offset);
986 else {
987 if (strncmp(origName, "FFT", 3) == 0)
988 offset = 3;
989 else if (strncmp(origName, "RealFFT", 7) == 0)
990 offset = 7;
991 else
992 offset = 0;
993 sprintf(s, "%s", origName + offset);
994 }
995 SDDS_CopyString(&name, s);
996 if (!origSymbol)
997 SDDS_CopyString(&origSymbol, origName + offset);
998 sprintf(s, "FFT %s", origSymbol);
999 SDDS_CopyString(&symbol, s);
1000
1001 sprintf(s, "Amplitude of FFT of %s", origSymbol);
1002 SDDS_CopyString(&description, s);
1003
1004 if (SDDS_NumberOfErrors() || (index0 = SDDS_DefineColumn(SDDSout, name, symbol, origUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
1005 return 0;
1006 free(name);
1007 free(symbol);
1008 free(description);
1009
1010 if (fftOffset == -1)
1011 fftOffset = 0;
1012
1013 if (psd_output & FL_PSDOUTPUT) {
1014 if (origUnits && !SDDS_StringIsBlank(origUnits)) {
1015 if (freqUnits && !SDDS_StringIsBlank(freqUnits)) {
1016 sprintf(s, "(%s)$a2$n/(%s)", origUnits, freqUnits);
1017 } else
1018 sprintf(s, "(%s)$a2$n", origUnits);
1019 SDDS_CopyString(&units, s);
1020 } else
1021 units = NULL;
1022
1023 sprintf(s, "PSD%s", origName + offset);
1024 SDDS_CopyString(&name, s);
1025
1026 if (!origSymbol)
1027 SDDS_CopyString(&origSymbol, origName + offset);
1028 sprintf(s, "PSD %s", origSymbol);
1029 SDDS_CopyString(&symbol, s);
1030
1031 sprintf(s, "PSD of %s", origSymbol);
1032 SDDS_CopyString(&description, s);
1033
1034 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, units, description, NULL, SDDS_DOUBLE, 0)) < 0)
1035 return 0;
1036 psdOffset = index1 - index0;
1037 free(name);
1038 if (units)
1039 free(units);
1040 free(symbol);
1041 free(description);
1042 }
1043
1044 if (psd_output & (FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT)) {
1045 if (origUnits && !SDDS_StringIsBlank(origUnits)) {
1046 SDDS_CopyString(&units, origUnits);
1047 } else
1048 units = NULL;
1049
1050 sprintf(s, "SqrtIntegPSD%s", origName + offset);
1051 SDDS_CopyString(&name, s);
1052
1053 if (!origSymbol)
1054 SDDS_CopyString(&origSymbol, origName + offset);
1055 sprintf(s, "Sqrt Integ PSD %s", origSymbol);
1056 SDDS_CopyString(&symbol, s);
1057
1058 sprintf(s, "Sqrt Integ PSD of %s", origSymbol);
1059 SDDS_CopyString(&description, s);
1060
1061 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, units, description, NULL, SDDS_DOUBLE, 0)) < 0)
1062 return 0;
1063 psdIntOffset = index1 - index0;
1064
1065 sprintf(s, "IntegPSD%s", origName + offset);
1066 SDDS_CopyString(&name, s);
1067
1068 if (!origSymbol)
1069 SDDS_CopyString(&origSymbol, origName + offset);
1070 sprintf(s, "Integ PSD %s", origSymbol);
1071 SDDS_CopyString(&symbol, s);
1072
1073 sprintf(s, "Integ PSD of %s", origSymbol);
1074 SDDS_CopyString(&description, s);
1075
1076 if (origUnits && !SDDS_StringIsBlank(origUnits)) {
1077 sprintf(s, "%sPower", origUnits);
1078 SDDS_CopyString(&units, origUnits);
1079 } else
1080 units = NULL;
1081
1082 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, units, description, NULL, SDDS_DOUBLE, 0)) < 0)
1083 return 0;
1084 psdIntPowerOffset = index1 - index0;
1085 free(name);
1086 if (units)
1087 free(units);
1088 free(symbol);
1089 free(description);
1090 }
1091
1092 if (full_output) {
1093 if (!inverse)
1094 sprintf(s, "RealFFT%s", origName + offset);
1095 else
1096 sprintf(s, "Real%s", origName + offset);
1097 SDDS_CopyString(&name, s);
1098
1099 if (!origSymbol)
1100 SDDS_CopyString(&origSymbol, origName + offset);
1101 if (!inverse)
1102 sprintf(s, "Re[FFT %s]", origSymbol);
1103 else
1104 sprintf(s, "Re[%s]", origSymbol);
1105 SDDS_CopyString(&symbol, s);
1106
1107 if (!inverse)
1108 sprintf(s, "Real part of FFT of %s", origSymbol);
1109 else
1110 sprintf(s, "Real part of %s", origSymbol);
1111 SDDS_CopyString(&description, s);
1112
1113 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, origUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
1114 return 0;
1115 realOffset = index1 - index0;
1116 free(name);
1117 free(symbol);
1118 free(description);
1119
1120 if (!inverse)
1121 sprintf(s, "ImagFFT%s", origName + offset);
1122 else
1123 sprintf(s, "Imag%s", origName + offset);
1124 SDDS_CopyString(&name, s);
1125
1126 if (!origSymbol)
1127 SDDS_CopyString(&origSymbol, origName + offset);
1128 if (!inverse)
1129 sprintf(s, "Im[FFT %s]", origSymbol);
1130 else
1131 sprintf(s, "Im[%s]", origSymbol);
1132 SDDS_CopyString(&symbol, s);
1133
1134 if (!inverse)
1135 sprintf(s, "Imaginary part of FFT of %s", origSymbol);
1136 else
1137 sprintf(s, "Imaginary part of %s", origSymbol);
1138 SDDS_CopyString(&description, s);
1139
1140 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, origUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
1141 return 0;
1142 imagOffset = index1 - index0;
1143 free(name);
1144 free(symbol);
1145 free(description);
1146
1147 if (!inverse)
1148 sprintf(s, "ArgFFT%s", origName + offset);
1149 else
1150 sprintf(s, "Arg%s", origName + offset);
1151 SDDS_CopyString(&name, s);
1152
1153 if (!origSymbol)
1154 SDDS_CopyString(&origSymbol, origName + offset);
1155 if (!inverse)
1156 sprintf(s, "Arg[FFT %s]", origSymbol);
1157 else
1158 sprintf(s, "Arg[%s]", origSymbol);
1159 SDDS_CopyString(&symbol, s);
1160
1161 if (!inverse)
1162 sprintf(s, "Phase of FFT of %s", origSymbol);
1163 else
1164 sprintf(s, "Phase of %s", origSymbol);
1165 SDDS_CopyString(&description, s);
1166
1167 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, "degrees", description, NULL, SDDS_DOUBLE, 0)) < 0)
1168 return 0;
1169 argOffset = index1 - index0;
1170 free(name);
1171 free(symbol);
1172 free(description);
1173 if (unwrap_phase) {
1174 if (!inverse)
1175 sprintf(s, "UnwrapArgFFT%s", origName + offset);
1176 else
1177 sprintf(s, "UnwrapArg%s", origName + offset);
1178 SDDS_CopyString(&name, s);
1179
1180 if (!origSymbol)
1181 SDDS_CopyString(&origSymbol, origName + offset);
1182 if (!inverse)
1183 sprintf(s, "UnwrapArg[FFT %s]", origSymbol);
1184 else
1185 sprintf(s, "UnwrapArg[%s]", origSymbol);
1186 SDDS_CopyString(&symbol, s);
1187
1188 if (!inverse)
1189 sprintf(s, "Unwrapped Phase of FFT of %s", origSymbol);
1190 else
1191 sprintf(s, "Unwrapped Phase of %s", origSymbol);
1192 SDDS_CopyString(&description, s);
1193
1194 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, "degrees", description, NULL, SDDS_DOUBLE, 0)) < 0)
1195 return 0;
1196 unwrappedArgOffset = index1 - index0;
1197 free(name);
1198 free(symbol);
1199 free(description);
1200 }
1201 }
1202
1203 return 1;
1204}
1205
1206void moveToStringArrayComplex(char ***targetReal, char ***targetImag, long *targets, char **sourceReal, char **sourceImag, long sources);
1207
1208long expandComplexColumnPairNames(SDDS_DATASET *SDDSin, char **name, char ***realName, char ***imagName, long names, char **excludeName, long excludeNames, long typeMode, long typeValue) {
1209 long i, j, k, realNames, imagNames, names2;
1210 char **realName1, **imagName1, **realName2, **imagName2;
1211 char *realPattern, *imagPattern = NULL;
1212 long longest=0;
1213
1214 if (!names || !name)
1215 return 0;
1216 realName1 = imagName1 = realName2 = imagName2 = NULL;
1217 realNames = imagNames = names2 = 0;
1218 for (i = 0; i < names; i++) {
1219 if (strlen(name[i]) > longest)
1220 longest = strlen(name[i]);
1221 }
1222 longest += 10;
1223 if (!(realPattern = SDDS_Malloc(sizeof(*realPattern) * longest)) || !(imagPattern = SDDS_Malloc(sizeof(*imagPattern) * longest)))
1224 SDDS_Bomb("memory allocation failure");
1225
1226 for (i = 0; i < names; i++) {
1227 for (j = 0; j < 2; j++) {
1228 if (j == 0) {
1229 sprintf(realPattern, "Real%s", name[i]);
1230 sprintf(imagPattern, "Imag%s", name[i]);
1231 } else {
1232 sprintf(realPattern, "%sReal", name[i]);
1233 sprintf(imagPattern, "%sImag", name[i]);
1234 }
1235 switch (typeMode) {
1236 case FIND_ANY_TYPE:
1237 case FIND_NUMERIC_TYPE:
1238 case FIND_INTEGER_TYPE:
1239 case FIND_FLOATING_TYPE:
1240 realNames = SDDS_MatchColumns(SDDSin, &realName1, SDDS_MATCH_STRING, typeMode, realPattern, SDDS_0_PREVIOUS | SDDS_OR);
1241 imagNames = SDDS_MatchColumns(SDDSin, &imagName1, SDDS_MATCH_STRING, typeMode, imagPattern, SDDS_0_PREVIOUS | SDDS_OR);
1242 break;
1243 case FIND_SPECIFIED_TYPE:
1244 if (!SDDS_VALID_TYPE(typeValue))
1245 SDDS_Bomb("invalid type value in expandColumnPairNames");
1246 realNames = SDDS_MatchColumns(SDDSin, &realName1, SDDS_MATCH_STRING, typeMode, typeValue, realPattern, SDDS_0_PREVIOUS | SDDS_OR);
1247 imagNames = SDDS_MatchColumns(SDDSin, &imagName1, SDDS_MATCH_STRING, typeMode, typeValue, imagPattern, SDDS_0_PREVIOUS | SDDS_OR);
1248 break;
1249 default:
1250 SDDS_Bomb("invalid typeMode in expandColumnPairNames");
1251 exit(EXIT_FAILURE);
1252 break;
1253 }
1254 if (realNames == 0)
1255 continue;
1256 if (realNames == -1 || imagNames == -1) {
1257 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
1258 SDDS_Bomb("unable to perform column name match in expandColumnPairNames");
1259 }
1260 if (realNames != imagNames)
1261 SDDS_Bomb("found different number of real and imaginary columns");
1262 if (excludeNames) {
1263 for (j = 0; j < excludeNames; j++)
1264 for (k = 0; k < realNames; k++)
1265 if (wild_match(realName1[k], excludeName[j])) {
1266 free(realName1[k]);
1267 free(imagName1[k]);
1268 imagName1[k] = realName1[k] = NULL;
1269 }
1270 }
1271 moveToStringArrayComplex(&realName2, &imagName2, &names2, realName1, imagName1, realNames);
1272 free(realName1);
1273 free(imagName1);
1274 }
1275 }
1276 free(realPattern);
1277 free(imagPattern);
1278 if (names2 == 0)
1279 return 0;
1280 *realName = realName2;
1281 *imagName = imagName2;
1282 return names2;
1283}
1284
1285void moveToStringArrayComplex(char ***targetReal, char ***targetImag, long *targets, char **sourceReal, char **sourceImag, long sources) {
1286 long i, j;
1287 if (!sources)
1288 return;
1289 if (!(*targetReal = SDDS_Realloc(*targetReal, sizeof(**targetReal) * (*targets + sources))) ||
1290 !(*targetImag = SDDS_Realloc(*targetImag, sizeof(**targetImag) * (*targets + sources))))
1291 SDDS_Bomb("memory allocation failure");
1292 for (i = 0; i < sources; i++) {
1293 if (sourceReal[i] == NULL || sourceImag[i] == NULL)
1294 continue;
1295 for (j = 0; j < *targets; j++)
1296 if (strcmp(sourceReal[i], (*targetReal)[j]) == 0)
1297 break;
1298 if (j == *targets) {
1299 (*targetReal)[j] = sourceReal[i];
1300 (*targetImag)[j] = sourceImag[i];
1301 *targets += 1;
1302 }
1303 }
1304}
SDDS.h
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
SDDS_CopyParameters
int32_t SDDS_CopyParameters(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
Definition SDDS_copy.c:286
SDDS_StartPage
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
Definition SDDS_dataprep.c:70
SDDS_SetParameters
int32_t SDDS_SetParameters(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Definition SDDS_dataprep.c:375
SDDS_SetColumn
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.
Definition SDDS_dataprep.c:1542
SDDS_AssertRowFlags
int32_t SDDS_AssertRowFlags(SDDS_DATASET *SDDS_dataset, uint32_t mode,...)
Sets acceptance flags for rows based on specified criteria.
Definition SDDS_extract.c:148
SDDS_GetParameterAsLong
int32_t * SDDS_GetParameterAsLong(SDDS_DATASET *SDDS_dataset, char *parameter_name, int32_t *memory)
Retrieves the value of a specified parameter as a 32-bit integer from the current data table of a dat...
Definition SDDS_extract.c:2616
SDDS_CountRowsOfInterest
int64_t SDDS_CountRowsOfInterest(SDDS_DATASET *SDDS_dataset)
Counts the number of rows marked as "of interest" in the current data table.
Definition SDDS_extract.c:364
SDDS_GetColumnInDoubles
double * SDDS_GetColumnInDoubles(SDDS_DATASET *SDDS_dataset, char *column_name)
Retrieves the data of a specified numerical column as an array of doubles, considering only rows mark...
Definition SDDS_extract.c:872
SDDS_GetColumnInformation
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.
Definition SDDS_info.c:41
SDDS_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:50
SDDS_Terminate
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1471
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1007
SDDS_InitializeOutput
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.
Definition SDDS_output.c:602
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1229
SDDS_DefineColumn
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.
Definition SDDS_output.c:1716
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:900
SDDS_DefineParameter
int32_t SDDS_DefineParameter(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, char *fixed_value)
Defines a data parameter with a fixed string value.
Definition SDDS_output.c:1473
SDDS_TransferAllParameterDefinitions
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.
Definition SDDS_transfer.c:343
SDDS_GetColumnIndex
int32_t SDDS_GetColumnIndex(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the index of a named column in the SDDS dataset.
Definition SDDS_utils.c:1375
SDDS_CheckColumn
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.
Definition SDDS_utils.c:4646
SDDS_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:474
SDDS_Malloc
void * SDDS_Malloc(size_t size)
Allocates memory of a specified size.
Definition SDDS_utils.c:705
SDDS_RegisterProgramName
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:318
SDDS_NumberOfErrors
int32_t SDDS_NumberOfErrors()
Retrieves the number of errors recorded by SDDS library routines.
Definition SDDS_utils.c:340
SDDS_StringIsBlank
int32_t SDDS_StringIsBlank(char *s)
Checks if a string is blank (contains only whitespace characters).
Definition SDDS_utils.c:2470
SDDS_MatchColumns
int32_t SDDS_MatchColumns(SDDS_DATASET *SDDS_dataset, char ***nameReturn, int32_t matchMode, int32_t typeMode,...)
Matches and retrieves column names from an SDDS dataset based on specified criteria.
Definition SDDS_utils.c:3559
SDDS_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:380
SDDS_CheckParameter
int32_t SDDS_CheckParameter(SDDS_DATASET *SDDS_dataset, char *name, char *units, int32_t type, FILE *fp_message)
Checks if a parameter exists in the SDDS dataset with the specified name, units, and type.
Definition SDDS_utils.c:4737
SDDS_CopyString
int32_t SDDS_CopyString(char **target, const char *source)
Copies a source string to a target string with memory allocation.
Definition SDDS_utils.c:922
SDDS_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:743
SDDS_VALID_TYPE
#define SDDS_VALID_TYPE(type)
Validates whether the given type identifier is within the defined range of SDDS types.
Definition SDDStypes.h:149
SDDS_STRING
#define SDDS_STRING
Identifier for the string data type.
Definition SDDStypes.h:85
SDDS_LONG
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
SDDS_ANY_NUMERIC_TYPE
#define SDDS_ANY_NUMERIC_TYPE
Special identifier used by SDDS_Check*() routines to accept any numeric type.
Definition SDDStypes.h:157
SDDS_DOUBLE
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
SDDSutils.h
Utility functions for SDDS dataset manipulation and string array operations.
tmalloc
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:65
largest_prime_factor
int64_t largest_prime_factor(int64_t number)
Find the largest prime factor of a number.
Definition factorize.c:83
find_min_max
int find_min_max(double *min, double *max, double *list, int64_t n)
Finds the minimum and maximum values in a list of doubles.
Definition findMinMax.c:33
ipow
double ipow(const double x, const int64_t p)
Compute x raised to the power p (x^p).
Definition ipow.c:33
match_string
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.
Definition match_string.c:34
compute_average
long compute_average(double *value, double *data, int64_t n)
Computes the average of an array of doubles.
Definition median.c:152
scanargs
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
Definition scanargs.c:36
processPipeOption
long processPipeOption(char **item, long items, unsigned long *flags)
Definition scanargs.c:357
processFilenames
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:391
scanItemList
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.
Definition scanitemlist.c:41
SDDS_DATASET
Definition SDDS.h:308
wild_match
int wild_match(char *string, char *template)
Determine whether one string is a wildcard match for another.
Definition wild_match.c:49
Generated by doxygen 1.12.0