SDDS ToolKit Programs and Libraries for C and Python
Loading...
Searching...
No Matches
sdds2dfft.c
Go to the documentation of this file.
1/**
2 * @file sdds2dfft.c
3 * @brief SDDS-format 2D FFT program.
4 *
5 * @details
6 * This program performs a two-dimensional Fast Fourier Transform (FFT) on data formatted in the
7 * Self Describing Data Set (SDDS) format. It provides various options for normalization, padding, truncation,
8 * suppressing averages, and more to customize the FFT process.
9 *
10 * @section Usage
11 * ```
12 * sdds2dfft [<inputfile>] [<outputfile>]
13 * [-pipe=[input][,output]]
14 * -columns=<indep-variable>[,<depen-quantity>[,...]]
15 * [-complexInput[=unfolded|folded]]
16 * [-exclude=<depen-quantity>[,...]]
17 * [-sampleInterval=<number>]
18 * [-normalize]
19 * [-fullOutput[=unfolded|folded],unwrapLimit=<value>]
20 * [-psdOutput[=plain][,{integrated|rintegrated[=<cutoff>]}]]
21 * [-inverse]
22 * [-padwithzeroes[=exponent]]
23 * [-truncate]
24 * [-suppressaverage]
25 * [-noWarnings]
26 * [-majorOrder=row|column]
27 * [-threads=<number>]
28 * ```
29 *
30 * @section Options
31 * | Required | Description |
32 * |---------------------------------------|---------------------------------------------------------------------------------------|
33 * | `-columns` | Specify independent and dependent variables for FFT analysis. |
34 *
35 * | Optional | Description |
36 * |---------------------|-----------------------------------------------------------------|
37 * | `-pipe` | Standard SDDS Toolkit pipe option. |
38 * | `-complexInput` | Indicates the input columns are in complex form. |
39 * | `-exclude` | Exclude quantities from analysis using wildcards. |
40 * | `-sampleInterval` | Request sampling of input data points at specified intervals. |
41 * | `-normalize` | Normalize output to a peak magnitude of 1. |
42 * | `-fullOutput` | Request real and imaginary parts of the FFT. |
43 * | `-psdOutput` | Request Power Spectral Density (PSD) output. |
44 * | `-inverse` | Perform inverse Fourier transform. |
45 * | `-padwithzeroes` | Pad data with zeroes to match required data points. |
46 * | `-truncate` | Truncate data to match required data points. |
47 * | `-suppressaverage` | Suppress the average value before FFT. |
48 * | `-noWarnings` | Suppress warning messages. |
49 * | `-majorOrder` | Specify output file's data order (row or column). |
50 * | `-threads` | Number of threads for row-wise FFT work. |
51 *
52 * @subsection Incompatibilities
53 * - `-inverse` is incompatible with:
54 * - `-complexInput=folded`
55 * - `-padwithzeroes` is incompatible with `-truncate`.
56 * - For `-complexInput`:
57 * - Requires `-columns` specifying dependent quantities in pairs (real, imaginary).
58 *
59 * @copyright
60 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
61 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
62 *
63 * @license
64 * This file is distributed under the terms of the Software License Agreement
65 * found in the file LICENSE included with this distribution.
66 *
67 * @authors
68 * H. Shang, R. Soliday
69 */
70
71#include "mdb.h"
72#include "SDDS.h"
73#include "scan.h"
74#include "fftpackC.h"
75#include "SDDSutils.h"
76#include <ctype.h>
77
78/* Enumeration for option types */
79enum option_type {
80 SET_NORMALIZE,
81 SET_PADWITHZEROES,
82 SET_TRUNCATE,
83 SET_SUPPRESSAVERAGE,
84 SET_SAMPLEINTERVAL,
85 SET_COLUMNS,
86 SET_FULLOUTPUT,
87 SET_PIPE,
88 SET_PSDOUTPUT,
89 SET_EXCLUDE,
90 SET_NOWARNINGS,
91 SET_COMPLEXINPUT,
92 SET_INVERSE,
93 SET_MAJOR_ORDER,
94 SET_THREADS,
95 N_OPTIONS
96};
97
98char *option[N_OPTIONS] = {
99 "normalize",
100 "padwithzeroes",
101 "truncate",
102 "suppressaverage",
103 "sampleinterval",
104 "columns",
105 "fulloutput",
106 "pipe",
107 "psdoutput",
108 "exclude",
109 "nowarnings",
110 "complexinput",
111 "inverse",
112 "majorOrder",
113 "threads",
114};
115
116#define FL_TRUNCATE 0x0001
117#define FL_PADWITHZEROES 0x0002
118#define FL_NORMALIZE 0x0004
119#define FL_SUPPRESSAVERAGE 0x0008
120#define FL_FULLOUTPUT 0x0010
121#define FL_MAKEFREQDATA 0x0020
122#define FL_PSDOUTPUT 0x0040
123#define FL_PSDINTEGOUTPUT 0x0080
124#define FL_PSDRINTEGOUTPUT 0x0100
125#define FL_FULLOUTPUT_FOLDED 0x0200
126#define FL_FULLOUTPUT_UNFOLDED 0x0400
127#define FL_COMPLEXINPUT_FOLDED 0x0800
128#define FL_COMPLEXINPUT_UNFOLDED 0x1000
129#define FL_UNWRAP_PHASE 0x2000
130
131static char *USAGE1 =
132 "Usage: sdds2dfft [<inputfile>] [<outputfile>]\n"
133 " [-pipe=[input][,output]]\n"
134 " -columns=<indep-variable>[,<depen-quantity>[,...]]\n"
135 " [-complexInput[=unfolded|folded]]\n"
136 " [-exclude=<depen-quantity>[,...]]\n"
137 " [-sampleInterval=<number>]\n"
138 " [-normalize]\n"
139 " [-fullOutput[=unfolded|folded],unwrapLimit=<value>]\n"
140 " [-psdOutput[=plain][,{integrated|rintegrated[=<cutoff>]}]]\n"
141 " [-inverse]\n"
142 " [-padwithzeroes[=exponent]]\n"
143 " [-truncate]\n"
144 " [-suppressaverage]\n"
145 " [-noWarnings]\n"
146 " [-majorOrder=row|column]\n"
147 " [-threads=<number>]\n"
148 "Options:\n"
149 " -pipe=[input][,output]\n"
150 " The standard SDDS Toolkit pipe option.\n"
151 " -columns=<indep-variable>[,<depen-quantity>[,...]]\n"
152 " Specifies the independent variable and dependent quantities to Fourier analyze.\n"
153 " <depen-quantity> entries may contain wildcards.\n"
154 " -complexInput[=unfolded|folded]\n"
155 " Indicates that the input columns are in complex form.\n"
156 " Options:\n"
157 " unfolded - The input frequency space is unfolded and must include negative frequencies.\n"
158 " folded - The input frequency space is folded (default).\n"
159 " -inverse\n"
160 " Produces the inverse Fourier transform. The output is always an unfolded spectrum.\n"
161 " If combined with -fullOutput=folded, it will be changed to -fullOutput=unfolded.\n"
162 " -exclude=<depen-quantity>[,...]\n"
163 " Specifies a list of wild-card patterns to exclude certain quantities from analysis.\n";
164
165static char *USAGE2 =
166 " -sampleInterval=<number>\n"
167 " Requests sampling of the input data points with the given interval.\n"
168 " -normalize\n"
169 " Normalizes the output to a peak magnitude of 1.\n"
170 " -fullOutput[=unfolded|folded],unwrapLimit=<value>\n"
171 " Requests output of the real and imaginary parts of the FFT.\n"
172 " Options:\n"
173 " unfolded - Outputs the unfolded frequency-space (full FFT).\n"
174 " folded - Outputs the folded frequency-space (half FFT) (default).\n"
175 " Additional parameter:\n"
176 " unwrapLimit=<value> - Unwraps the phase where the relative magnitude exceeds this limit.\n"
177 " -psdOutput[=plain][,{integrated|rintegrated[=<cutoff>]}]\n"
178 " Requests output of the Power Spectral Density (PSD).\n"
179 " Qualifiers:\n"
180 " plain - Includes plain PSD output.\n"
181 " integrated - Includes integrated PSD.\n"
182 " rintegrated - Includes reverse-integrated PSD with an optional cutoff frequency.\n"
183 " -padwithzeroes[=exponent] | -truncate\n"
184 " -padwithzeroes: Pads the data with zeroes if the number of data points is not a product of small primes.\n"
185 " Optionally specify an exponent to determine the padding factor.\n"
186 " -truncate: Truncates the data if the number of data points is not a product of small primes.\n"
187 " -suppressaverage\n"
188 " Removes the average value of the data before performing the FFT.\n"
189 " -noWarnings\n"
190 " Suppresses warning messages.\n"
191 " -majorOrder=row|column\n"
192 " Specifies the output file's data order.\n"
193 " row - Row-major order.\n"
194 " column - Column-major order.\n"
195 " -threads=<number>\n"
196 " Number of threads for row-wise FFT work.\n\n"
197 "Program by Hairong Shang. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n";
198
199int64_t greatestProductOfSmallPrimes(int64_t rows);
200long create_fft_frequency_column(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *timeName, char *freqUnits, long inverse);
201
202long create_fft_columns(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *origName, char *indepName,
203 char *freqUnits, long full_output, unsigned long psd_output, long complexInput,
204 long inverse, long unwrap_phase);
205long create_fft_parameters(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *indepName, char *freqUnits);
206char *makeFrequencyUnits(SDDS_DATASET *SDDSin, char *indepName);
207long expandComplexColumnPairNames(SDDS_DATASET *SDDSin, char **name, char ***realName, char ***imagName,
208 long names, char **excludeName, long excludeNames, long typeMode, long typeValue);
209void moveToStringArrayComplex(char ***targetReal, char ***targetImag, long *targets, char **sourceReal, char **sourceImag, long sources);
210
211int main(int argc, char **argv) {
212 int iArg, j;
213 char *freqUnits;
214 char *indepQuantity, **depenQuantity, **exclude, **realQuan = NULL, **imagQuan = NULL;
215 long depenQuantities, excludes;
216 char *input, *output;
217 long sampleInterval, readCode, noWarnings, complexInput, inverse, spectrumFoldParExist = 0, colsToUse;
218 int64_t i, rows, rowsToUse, primeRows, pow2Rows, n_freq, fftrows;
219 int32_t spectrumFolded = 0, page = 0, index;
220 unsigned long flags, pipeFlags, complexInputFlags = 0, fullOutputFlags = 0, majorOrderFlag;
221 long primeCols, pow2Cols;
222 SCANNED_ARG *scanned;
223 SDDS_DATASET SDDSin, SDDSout;
224 double *tdata, rintegCutOffFreq, unwrapLimit = 0;
225 long padFactor;
226 short columnMajorOrder = -1;
227 int threads = 1;
228 double length, *real_imag = NULL, **real = NULL, **imag = NULL, *real_imag1 = NULL, *fdata = NULL, df, t0, factor;
229 double dtf_real, dtf_imag, *arg = NULL, *magData = NULL;
230 char str[256], *tempStr = NULL;
231
232 SDDS_RegisterProgramName(argv[0]);
233 argc = scanargs(&scanned, argc, argv);
234 if (argc < 3 || argc > (3 + N_OPTIONS)) {
235 fprintf(stderr, "%s%s", USAGE1, USAGE2);
236 exit(EXIT_FAILURE);
237 }
238 rintegCutOffFreq = 0;
239 output = input = NULL;
240 flags = pipeFlags = excludes = complexInput = inverse = 0;
241 sampleInterval = 1;
242 indepQuantity = NULL;
243 depenQuantity = exclude = NULL;
244 depenQuantities = 0;
245 noWarnings = 0;
246 padFactor = 0;
247
248 for (iArg = 1; iArg < argc; iArg++) {
249 if (scanned[iArg].arg_type == OPTION) {
250 /* process options here */
251 switch (match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
252 case SET_NORMALIZE:
253 flags |= FL_NORMALIZE;
254 break;
255 case SET_PADWITHZEROES:
256 flags |= FL_PADWITHZEROES;
257 if (scanned[iArg].n_items != 1) {
258 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%ld", &padFactor) != 1 || padFactor < 1)
259 SDDS_Bomb("invalid -padwithzeroes syntax");
260 }
261 break;
262 case SET_TRUNCATE:
263 flags |= FL_TRUNCATE;
264 break;
265 case SET_SUPPRESSAVERAGE:
266 flags |= FL_SUPPRESSAVERAGE;
267 break;
268 case SET_SAMPLEINTERVAL:
269 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%ld", &sampleInterval) != 1 || sampleInterval <= 0)
270 SDDS_Bomb("invalid -sampleinterval syntax");
271 break;
272 case SET_COLUMNS:
273 if (indepQuantity)
274 SDDS_Bomb("only one -columns option may be given");
275 if (scanned[iArg].n_items < 2)
276 SDDS_Bomb("invalid -columns syntax");
277 indepQuantity = scanned[iArg].list[1];
278 if (scanned[iArg].n_items >= 2) {
279 depenQuantity = tmalloc(sizeof(*depenQuantity) * (depenQuantities = scanned[iArg].n_items - 2));
280 for (i = 0; i < depenQuantities; i++)
281 SDDS_CopyString(&depenQuantity[i], scanned[iArg].list[i + 2]);
282 }
283 break;
284 case SET_FULLOUTPUT:
285 flags |= FL_FULLOUTPUT;
286 if (scanned[iArg].n_items >= 2) {
287 scanned[iArg].n_items--;
288 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))
289 SDDS_Bomb("Invalid -fullOutput syntax");
290 scanned[iArg].n_items++;
291 if (fullOutputFlags & FL_FULLOUTPUT_UNFOLDED)
292 flags |= FL_FULLOUTPUT_UNFOLDED;
293 else
294 flags |= FL_FULLOUTPUT_FOLDED;
295 if (fullOutputFlags & FL_UNWRAP_PHASE)
296 flags |= FL_UNWRAP_PHASE;
297 }
298 break;
299 case SET_PIPE:
300 if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
301 SDDS_Bomb("invalid -pipe syntax");
302 break;
303 case SET_PSDOUTPUT:
304 if (scanned[iArg].n_items > 1) {
305 unsigned long tmpFlags;
306 if (strchr(scanned[iArg].list[1], '=') == NULL) {
307 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))
308 SDDS_Bomb("invalid -psdOutput syntax");
309 } else {
310 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))
311 SDDS_Bomb("invalid -psdOutput syntax");
312 }
313 flags |= tmpFlags;
314 } else {
315 flags |= FL_PSDOUTPUT;
316 }
317 if ((flags & FL_PSDINTEGOUTPUT) && (flags & FL_PSDRINTEGOUTPUT))
318 SDDS_Bomb("invalid -psdOutput syntax: give only one of integrated or rintegrated");
319 break;
320 case SET_EXCLUDE:
321 if (scanned[iArg].n_items < 2)
322 SDDS_Bomb("invalid -exclude syntax");
323 for (j = 1; j < scanned[iArg].n_items; j++)
324 excludes = appendToStringArray(&exclude, excludes, scanned[iArg].list[j]);
325 break;
326 case SET_NOWARNINGS:
327 noWarnings = 1;
328 break;
329 case SET_COMPLEXINPUT:
330 complexInput = 1;
331 if (scanned[iArg].n_items == 2) {
332 scanned[iArg].n_items--;
333 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))
334 SDDS_Bomb("Invalid -complexInput syntax");
335 scanned[iArg].n_items++;
336 }
337 break;
338 case SET_INVERSE:
339 inverse = 1;
340 break;
341 case SET_MAJOR_ORDER:
342 majorOrderFlag = 0;
343 scanned[iArg].n_items--;
344 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)))
345 SDDS_Bomb("invalid -majorOrder syntax/values");
346 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
347 columnMajorOrder = 1;
348 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
349 columnMajorOrder = 0;
350 break;
351 case SET_THREADS:
352 if (scanned[iArg].n_items != 2 ||
353 sscanf(scanned[iArg].list[1], "%d", &threads) != 1 || threads < 1)
354 SDDS_Bomb("invalid -threads syntax");
355 break;
356 default:
357 fprintf(stderr, "error: unknown/ambiguous option: %s\n", scanned[iArg].list[0]);
358 exit(EXIT_FAILURE);
359 break;
360 }
361 } else {
362 if (!input)
363 input = scanned[iArg].list[0];
364 else if (!output)
365 output = scanned[iArg].list[0];
366 else
367 SDDS_Bomb("too many filenames seen");
368 }
369 }
370 if (!complexInput) {
371 if (!noWarnings && inverse)
372 fprintf(stderr, "Warning: The inverse option is ignored since it only works with -complexInput.\n");
373 inverse = 0;
374 }
375 if (!noWarnings && inverse && (flags & FL_FULLOUTPUT_FOLDED))
376 fprintf(stderr, "Warning: The combination of -inverse and -fullOutput=folded will be changed to -inverse -fullOutput=unfolded.\n");
377
378 processFilenames("sdds2dfft", &input, &output, pipeFlags, 0, NULL);
379
380 if (!indepQuantity)
381 SDDS_Bomb("Supply the independent quantity name with the -columns option");
382
383 if ((flags & FL_TRUNCATE) && (flags & FL_PADWITHZEROES))
384 SDDS_Bomb("Specify only one of -padwithzeroes and -truncate");
385 if (!inverse) {
386 /* For 2D FFT, always use full output unfolded */
387 flags |= FL_FULLOUTPUT;
388 flags |= FL_FULLOUTPUT_UNFOLDED;
389 }
390 if (!SDDS_InitializeInput(&SDDSin, input))
391 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
392
393 if (SDDS_CheckColumn(&SDDSin, indepQuantity, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OKAY)
394 exit(EXIT_FAILURE);
395
396 excludes = appendToStringArray(&exclude, excludes, indepQuantity);
397 if (!depenQuantities)
398 depenQuantities = appendToStringArray(&depenQuantity, depenQuantities, "*");
399
400 if (!complexInput) {
401 if ((depenQuantities = expandColumnPairNames(&SDDSin, &depenQuantity, NULL, depenQuantities, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
402 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
403 SDDS_Bomb("No quantities selected to FFT");
404 }
405 } else {
406 if ((depenQuantities = expandComplexColumnPairNames(&SDDSin, depenQuantity, &realQuan, &imagQuan, depenQuantities, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
407 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
408 SDDS_Bomb("No quantities selected to FFT");
409 }
410 }
411
412#if 0
413 fprintf(stderr, "%ld dependent quantities:\n", depenQuantities);
414 for (i = 0; i < depenQuantities; i++)
415 fprintf(stderr, " %s\n", depenQuantity[i]);
416#endif
417
418 if (!(freqUnits = makeFrequencyUnits(&SDDSin, indepQuantity)) ||
419 !SDDS_InitializeOutput(&SDDSout, SDDS_BINARY, 0, NULL, "sdds2dfft output", output) ||
420 !create_fft_frequency_column(&SDDSout, &SDDSin, indepQuantity, freqUnits, inverse) ||
421 SDDS_DefineParameter(&SDDSout, "fftFrequencies", NULL, NULL, NULL, NULL, SDDS_LONG, NULL) < 0 ||
422 SDDS_DefineParameter(&SDDSout, "fftFrequencySpacing", "$gD$rf", freqUnits, NULL, NULL, SDDS_DOUBLE, NULL) < 0)
423 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
424 if (columnMajorOrder != -1)
425 SDDSout.layout.data_mode.column_major = columnMajorOrder;
426 else
427 SDDSout.layout.data_mode.column_major = SDDSin.layout.data_mode.column_major;
428
429 if ((flags & FL_FULLOUTPUT) && SDDS_DefineParameter(&SDDSout, "SpectrumFolded", NULL, NULL, NULL, NULL, SDDS_LONG, NULL) < 0)
430 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
431 if (complexInput) {
432 if (!complexInputFlags) {
433 if (SDDS_CheckParameter(&SDDSin, "SpectrumFolded", NULL, SDDS_LONG, NULL) == SDDS_CHECK_OK)
434 spectrumFoldParExist = 1;
435 } else if (complexInputFlags & FL_COMPLEXINPUT_UNFOLDED)
436 flags |= FL_COMPLEXINPUT_UNFOLDED;
437 else
438 flags |= FL_COMPLEXINPUT_FOLDED;
439 }
440 for (i = 0; i < depenQuantities; i++) {
441 if (!complexInput)
442 create_fft_columns(&SDDSout, &SDDSin, depenQuantity[i], indepQuantity, freqUnits,
443 flags & FL_FULLOUTPUT, flags & (FL_PSDOUTPUT + FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT),
444 0, inverse, flags & FL_UNWRAP_PHASE);
445 else
446 create_fft_columns(&SDDSout, &SDDSin, realQuan[i], indepQuantity, freqUnits,
447 flags & FL_FULLOUTPUT, flags & (FL_PSDOUTPUT + FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT),
448 1, inverse, flags & FL_UNWRAP_PHASE);
449 }
450
451 if (!SDDS_TransferAllParameterDefinitions(&SDDSout, &SDDSin, SDDS_TRANSFER_KEEPOLD) || !SDDS_WriteLayout(&SDDSout))
452 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
453
454 colsToUse = depenQuantities;
455 primeCols = greatestProductOfSmallPrimes(depenQuantities);
456 if (depenQuantities != primeCols || padFactor) {
457 if (flags & FL_PADWITHZEROES) {
458 pow2Cols = ipow(2., ((long)(log((double)depenQuantities) / log(2.0F))) + (padFactor ? padFactor : 1));
459 if ((primeCols = greatestProductOfSmallPrimes(pow2Cols)) > depenQuantities)
460 colsToUse = primeCols;
461 else
462 colsToUse = pow2Cols;
463 fprintf(stdout, "Using %ld columns\n", colsToUse);
464 } else if (flags & FL_TRUNCATE)
465 colsToUse = greatestProductOfSmallPrimes(depenQuantities);
466 else if (largest_prime_factor(depenQuantities) > 100 && !noWarnings)
467 fputs("Warning: Number of dependent columns has large prime factors.\nThis could take a very long time.\nConsider using the -truncate option.\n", stderr);
468 }
469 real_imag = tmalloc(sizeof(*real_imag) * (2 * colsToUse + 2));
470 real = malloc(sizeof(*real) * colsToUse);
471 imag = malloc(sizeof(*imag) * colsToUse);
472
473 while ((readCode = SDDS_ReadPage(&SDDSin)) > 0) {
474 page++;
475 if ((rows = SDDS_CountRowsOfInterest(&SDDSin)) < 0)
476 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
477 if (page == 1 && spectrumFoldParExist) {
478 if (!SDDS_GetParameterAsLong(&SDDSin, "SpectrumFolded", &spectrumFolded))
479 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
480 if (spectrumFolded)
481 flags |= FL_COMPLEXINPUT_FOLDED;
482 else
483 flags |= FL_COMPLEXINPUT_UNFOLDED;
484 }
485 if (rows) {
486 rowsToUse = rows;
487 primeRows = greatestProductOfSmallPrimes(rows);
488 if (rows != primeRows || padFactor) {
489 if (flags & FL_PADWITHZEROES) {
490 pow2Rows = ipow(2., ((long)(log((double)rows) / log(2.0F))) + (padFactor ? padFactor : 1));
491 if ((primeRows = greatestProductOfSmallPrimes(pow2Rows)) > rows)
492 rowsToUse = primeRows;
493 else
494 rowsToUse = pow2Rows;
495 fprintf(stdout, "Using %" PRId64 " rows\n", rowsToUse);
496 } else if (flags & FL_TRUNCATE)
497 rowsToUse = greatestProductOfSmallPrimes(rows);
498 else if (largest_prime_factor(rows) > 100 && !noWarnings)
499 fputs("Warning: Number of points has large prime factors.\nThis could take a very long time.\nConsider using the -truncate option.\n", stderr);
500 }
501 if (!(tdata = SDDS_GetColumnInDoubles(&SDDSin, indepQuantity)))
502 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
503
504 for (j = 0; j < colsToUse; j++) {
505 real[j] = imag[j] = NULL;
506 if (j < depenQuantities) {
507 if (complexInput) {
508 if (!(real[j] = (double *)SDDS_GetColumnInDoubles(&SDDSin, realQuan[j])) ||
509 !(imag[j] = (double *)SDDS_GetColumnInDoubles(&SDDSin, imagQuan[j])))
510 SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
511 } else {
512 if (!(real[j] = (double *)SDDS_GetColumnInDoubles(&SDDSin, depenQuantity[j])))
513 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
514 imag[j] = calloc(sizeof(**imag), rowsToUse);
515 }
516 if (rows < rowsToUse) {
517 real[j] = SDDS_Realloc(real[j], sizeof(**real) * rowsToUse);
518 imag[j] = SDDS_Realloc(imag[j], sizeof(**imag) * rowsToUse);
519 }
520 } else {
521 real[j] = calloc(sizeof(**real), rowsToUse);
522 imag[j] = calloc(sizeof(**imag), rowsToUse);
523 }
524 }
525 fdata = malloc(sizeof(*fdata) * rowsToUse);
526 if (rows < rowsToUse) {
527 length = ((double)rows) * (tdata[rows - 1] - tdata[0]) / ((double)rows - 1.0);
528 } else
529 length = tdata[rows - 1] - tdata[0];
530 t0 = tdata[0];
531 df = factor = 1.0 / length;
532 free(tdata);
533 for (i = 0; i < rows; i++)
534 fdata[i] = i * df;
535 for (i = rows; i < rowsToUse; i++) {
536 fdata[i] = i * df;
537 }
538 /* First perform FFT per row */
539 if (threads <= 1) {
540 for (i = 0; i < rows; i++) {
541 for (j = 0; j < colsToUse; j++) {
542 real_imag[2 * j] = real_imag[2 * j + 1] = 0;
543 if (j < depenQuantities) {
544 real_imag[2 * j] = real[j][i];
545 if (imag[j])
546 real_imag[2 * j + 1] = imag[j][i];
547 else
548 real_imag[2 * j + 1] = 0;
549 }
550 }
551 complexFFT(real_imag, colsToUse, inverse);
552 for (j = 0; j < colsToUse; j++) {
553 real[j][i] = real_imag[2 * j];
554 imag[j][i] = real_imag[2 * j + 1];
555 }
556 }
557 } else {
558#pragma omp parallel if (threads > 1) num_threads(threads)
559 {
560 double *real_imag_thread = tmalloc(sizeof(*real_imag_thread) * (2 * colsToUse + 2));
561 int64_t row;
562#pragma omp for private(j)
563 for (row = 0; row < rows; row++) {
564 for (j = 0; j < colsToUse; j++) {
565 real_imag_thread[2 * j] = real_imag_thread[2 * j + 1] = 0;
566 if (j < depenQuantities) {
567 real_imag_thread[2 * j] = real[j][row];
568 if (imag[j])
569 real_imag_thread[2 * j + 1] = imag[j][row];
570 else
571 real_imag_thread[2 * j + 1] = 0;
572 }
573 }
574 complexFFT(real_imag_thread, colsToUse, inverse);
575 for (j = 0; j < colsToUse; j++) {
576 real[j][row] = real_imag_thread[2 * j];
577 imag[j][row] = real_imag_thread[2 * j + 1];
578 }
579 }
580 free(real_imag_thread);
581 }
582 }
583 /* Then perform FFT by column */
584 n_freq = rowsToUse;
585 fftrows = rowsToUse;
586 arg = malloc(sizeof(*arg) * rowsToUse);
587 magData = malloc(sizeof(*magData) * rowsToUse);
588 real_imag1 = calloc(sizeof(*real_imag1), 2 * fftrows + 2);
589
590 for (j = 0; j < depenQuantities; j++) {
591 for (i = 0; i < rowsToUse; i++) {
592 if (i < rows) {
593 real_imag1[2 * i] = real[j][i];
594 real_imag1[2 * i + 1] = imag[j][i];
595 } else {
596 real_imag1[2 * i] = 0;
597 real_imag1[2 * i + 1] = 0;
598 }
599 }
600 complexFFT(real_imag1, rowsToUse, inverse);
601 for (i = 0; i < n_freq; i++) {
602 dtf_real = cos(-2 * PI * fdata[i] * t0);
603 dtf_imag = sin(-2 * PI * fdata[i] * t0);
604 real[j][i] = real_imag1[2 * i] * dtf_real - real_imag1[2 * i + 1] * dtf_imag;
605 imag[j][i] = real_imag1[2 * i + 1] * dtf_real + real_imag1[2 * i] * dtf_imag;
606 magData[i] = sqrt(sqr(real[j][i]) + sqr(imag[j][i]));
607 if (real[j][i] || imag[j][i])
608 arg[i] = 180.0 / PI * atan2(imag[j][i], real[j][i]);
609 else
610 arg[i] = 0;
611 }
612 if (flags & FL_NORMALIZE) {
613 factor = -DBL_MAX;
614 for (i = 0; i < n_freq; i++)
615 if (magData[i] > factor)
616 factor = magData[i];
617 if (factor != -DBL_MAX)
618 for (i = 0; i < n_freq; i++) {
619 real[j][i] /= factor;
620 imag[j][i] /= factor;
621 magData[i] /= factor;
622 }
623 }
624 if (!inverse)
625 sprintf(str, "FFT%s", depenQuantity[j] + (imagQuan ? 4 : 0));
626 else {
627 if (complexInput)
628 tempStr = realQuan[j];
629 else
630 tempStr = depenQuantity[j];
631
632 if (strncmp(tempStr, "FFT", 3) == 0)
633 sprintf(str, "%s", tempStr + 3);
634 else if (strncmp(tempStr, "RealFFT", 7) == 0)
635 sprintf(str, "%s", tempStr + 7);
636 else
637 sprintf(str, "%s", tempStr);
638 }
639 if ((index = SDDS_GetColumnIndex(&SDDSout, str)) < 0)
640 exit(EXIT_FAILURE);
641 if (!SDDS_StartPage(&SDDSout, rowsToUse) ||
642 !SDDS_CopyParameters(&SDDSout, &SDDSin) ||
643 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, "fftFrequencies", n_freq, "fftFrequencySpacing", df, NULL))
644 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
645 if (flags & FL_FULLOUTPUT) {
646 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, magData, n_freq, index) ||
647 !SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, real[j], n_freq, index + 1) ||
648 !SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, imag[j], n_freq, index + 2) ||
649 !SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, arg, n_freq, index + 3))
650 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
651 } else {
652 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, real[j], n_freq, index))
653 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
654 }
655 }
656 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, fdata, n_freq, 0))
657 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
658 free(fdata);
659 free(arg);
660 free(magData);
661 for (j = 0; j < colsToUse; j++) {
662 if (real[j])
663 free(real[j]);
664 if (imag[j])
665 free(imag[j]);
666 }
667 free(real_imag1);
668 } else {
669 if (!SDDS_StartPage(&SDDSout, 0) || !SDDS_CopyParameters(&SDDSout, &SDDSin))
670 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
671 }
672 if (!SDDS_WritePage(&SDDSout))
673 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
674 }
675 if (!SDDS_Terminate(&SDDSin)) {
676 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
677 exit(EXIT_FAILURE);
678 }
679 if (!SDDS_Terminate(&SDDSout)) {
680 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
681 exit(EXIT_FAILURE);
682 }
683 if (excludes) {
684 SDDS_FreeStringArray(exclude, excludes);
685 free(exclude);
686 }
687 free(real);
688 free(imag);
689 if (realQuan) {
690 SDDS_FreeStringArray(realQuan, depenQuantities);
691 SDDS_FreeStringArray(imagQuan, depenQuantities);
692 free(realQuan);
693 free(imagQuan);
694 } else {
695 SDDS_FreeStringArray(depenQuantity, depenQuantities);
696 free(depenQuantity);
697 }
698 free(real_imag);
699 free_scanargs(&scanned, argc);
700 return EXIT_SUCCESS;
701}
702
703static long psdOffset, argOffset, realOffset, imagOffset, fftOffset = -1, psdIntOffset, unwrappedArgOffset = -1;
704
705long create_fft_frequency_column(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *timeName, char *freqUnits, long inverse) {
706 char s[SDDS_MAXLINE];
707 char *timeSymbol;
708 char *description;
709
710 if (SDDS_GetColumnInformation(SDDSin, "symbol", &timeSymbol, SDDS_GET_BY_NAME, timeName) != SDDS_STRING)
711 return 0;
712 if (!timeSymbol || SDDS_StringIsBlank(timeSymbol))
713 SDDS_CopyString(&timeSymbol, timeName);
714
715 sprintf(s, "Frequency for %s", timeSymbol);
716 SDDS_CopyString(&description, s);
717 if (!inverse) {
718 if (SDDS_DefineColumn(SDDSout, "f", NULL, freqUnits, description, NULL, SDDS_DOUBLE, 0) < 0) {
719 free(timeSymbol);
720 free(description);
721 return 0;
722 }
723 } else {
724 sprintf(s, "inverse for %s", timeSymbol);
725 SDDS_CopyString(&description, s);
726 if (SDDS_DefineColumn(SDDSout, "t", NULL, freqUnits, description, NULL, SDDS_DOUBLE, 0) < 0) {
727 free(timeSymbol);
728 free(description);
729 return 0;
730 }
731 }
732 free(timeSymbol);
733 free(description);
734 return 1;
735}
736
737long create_fft_columns(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *origName, char *indepName,
738 char *freqUnits, long full_output, unsigned long psd_output, long complexInput,
739 long inverse, long unwrap_phase) {
740 char s[SDDS_MAXLINE];
741 char *origUnits, *origSymbol;
742 char *description, *name, *symbol, *units;
743 long index0, index1;
744 long offset = 0;
745
746 if (complexInput)
747 offset = 4;
748 if (SDDS_GetColumnInformation(SDDSin, "units", &origUnits, SDDS_GET_BY_NAME, origName) != SDDS_STRING ||
749 SDDS_GetColumnInformation(SDDSin, "symbol", &origSymbol, SDDS_GET_BY_NAME, origName) != SDDS_STRING)
750 return 0;
751 if (!inverse)
752 sprintf(s, "FFT%s", origName + offset);
753 else {
754 if (strncmp(origName, "FFT", 3) == 0)
755 offset = 3;
756 else if (strncmp(origName, "RealFFT", 7) == 0)
757 offset = 7;
758 else
759 offset = 0;
760 sprintf(s, "%s", origName + offset);
761 }
762 SDDS_CopyString(&name, s);
763 if (!origSymbol)
764 SDDS_CopyString(&origSymbol, origName + offset);
765 sprintf(s, "FFT %s", origSymbol);
766 SDDS_CopyString(&symbol, s);
767
768 sprintf(s, "Amplitude of FFT of %s", origSymbol);
769 SDDS_CopyString(&description, s);
770
771 if (SDDS_NumberOfErrors() ||
772 (index0 = SDDS_DefineColumn(SDDSout, name, symbol, origUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
773 return 0;
774 free(name);
775 free(symbol);
776 free(description);
777
778 if (fftOffset == -1)
779 fftOffset = 0;
780
781 if (psd_output & FL_PSDOUTPUT) {
782 if (origUnits && !SDDS_StringIsBlank(origUnits)) {
783 if (freqUnits && !SDDS_StringIsBlank(freqUnits)) {
784 sprintf(s, "(%s)$a2$n/(%s)", origUnits, freqUnits);
785 } else
786 sprintf(s, "(%s)$a2$n", origUnits);
787 SDDS_CopyString(&units, s);
788 } else
789 units = NULL;
790
791 sprintf(s, "PSD%s", origName + offset);
792 SDDS_CopyString(&name, s);
793
794 if (!origSymbol)
795 SDDS_CopyString(&origSymbol, origName + offset);
796 sprintf(s, "PSD %s", origSymbol);
797 SDDS_CopyString(&symbol, s);
798
799 sprintf(s, "PSD of %s", origSymbol);
800 SDDS_CopyString(&description, s);
801
802 if (SDDS_NumberOfErrors() ||
803 (index1 = SDDS_DefineColumn(SDDSout, name, symbol, units, description, NULL, SDDS_DOUBLE, 0)) < 0)
804 return 0;
805 psdOffset = index1 - index0;
806 free(name);
807 if (units)
808 free(units);
809 free(symbol);
810 free(description);
811 }
812
813 if (psd_output & (FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT)) {
814 if (origUnits && !SDDS_StringIsBlank(origUnits)) {
815 SDDS_CopyString(&units, origUnits);
816 } else
817 units = NULL;
818
819 sprintf(s, "SqrtIntegPSD%s", origName + offset);
820 SDDS_CopyString(&name, s);
821
822 if (!origSymbol)
823 SDDS_CopyString(&origSymbol, origName + offset);
824 sprintf(s, "Sqrt Integ PSD %s", origSymbol);
825 SDDS_CopyString(&symbol, s);
826
827 sprintf(s, "Sqrt Integ PSD of %s", origSymbol);
828 SDDS_CopyString(&description, s);
829
830 if (SDDS_NumberOfErrors() ||
831 (index1 = SDDS_DefineColumn(SDDSout, name, symbol, units, description, NULL, SDDS_DOUBLE, 0)) < 0)
832 return 0;
833 psdIntOffset = index1 - index0;
834 free(name);
835 if (units)
836 free(units);
837 free(symbol);
838 free(description);
839 }
840
841 if (full_output) {
842 if (!inverse)
843 sprintf(s, "RealFFT%s", origName + offset);
844 else
845 sprintf(s, "Real%s", origName + offset);
846 SDDS_CopyString(&name, s);
847
848 if (!origSymbol)
849 SDDS_CopyString(&origSymbol, origName + offset);
850 if (!inverse)
851 sprintf(s, "Re[FFT %s]", origSymbol);
852 else
853 sprintf(s, "Re[%s]", origSymbol);
854 SDDS_CopyString(&symbol, s);
855
856 if (!inverse)
857 sprintf(s, "Real part of FFT of %s", origSymbol);
858 else
859 sprintf(s, "Real part of %s", origSymbol);
860 SDDS_CopyString(&description, s);
861
862 if (SDDS_NumberOfErrors() ||
863 (index1 = SDDS_DefineColumn(SDDSout, name, symbol, origUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
864 return 0;
865 realOffset = index1 - index0;
866 free(name);
867 free(symbol);
868 free(description);
869
870 if (!inverse)
871 sprintf(s, "ImagFFT%s", origName + offset);
872 else
873 sprintf(s, "Imag%s", origName + offset);
874 SDDS_CopyString(&name, s);
875
876 if (!origSymbol)
877 SDDS_CopyString(&origSymbol, origName + offset);
878 if (!inverse)
879 sprintf(s, "Im[FFT %s]", origSymbol);
880 else
881 sprintf(s, "Im[%s]", origSymbol);
882 SDDS_CopyString(&symbol, s);
883
884 if (!inverse)
885 sprintf(s, "Imaginary part of FFT of %s", origSymbol);
886 else
887 sprintf(s, "Imaginary part of %s", origSymbol);
888 SDDS_CopyString(&description, s);
889
890 if (SDDS_NumberOfErrors() ||
891 (index1 = SDDS_DefineColumn(SDDSout, name, symbol, origUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
892 return 0;
893 imagOffset = index1 - index0;
894 free(name);
895 free(symbol);
896 free(description);
897
898 if (!inverse)
899 sprintf(s, "ArgFFT%s", origName + offset);
900 else
901 sprintf(s, "Arg%s", origName + offset);
902 SDDS_CopyString(&name, s);
903
904 if (!origSymbol)
905 SDDS_CopyString(&origSymbol, origName + offset);
906 if (!inverse)
907 sprintf(s, "Arg[FFT %s]", origSymbol);
908 else
909 sprintf(s, "Arg[%s]", origSymbol);
910 SDDS_CopyString(&symbol, s);
911
912 if (!inverse)
913 sprintf(s, "Phase of FFT of %s", origSymbol);
914 else
915 sprintf(s, "Phase of %s", origSymbol);
916 SDDS_CopyString(&description, s);
917
918 if (SDDS_NumberOfErrors() ||
919 (index1 = SDDS_DefineColumn(SDDSout, name, symbol, "degrees", description, NULL, SDDS_DOUBLE, 0)) < 0)
920 return 0;
921 argOffset = index1 - index0;
922 free(name);
923 free(symbol);
924 free(description);
925 if (unwrap_phase) {
926 if (!inverse)
927 sprintf(s, "UnwrapArgFFT%s", origName + offset);
928 else
929 sprintf(s, "UnwrapArg%s", origName + offset);
930 SDDS_CopyString(&name, s);
931
932 if (!origSymbol)
933 SDDS_CopyString(&origSymbol, origName + offset);
934 if (!inverse)
935 sprintf(s, "UnwrapArg[FFT %s]", origSymbol);
936 else
937 sprintf(s, "UnwrapArg[%s]", origSymbol);
938 SDDS_CopyString(&symbol, s);
939
940 if (!inverse)
941 sprintf(s, "Unwrapped Phase of FFT of %s", origSymbol);
942 else
943 sprintf(s, "Unwrapped Phase of %s", origSymbol);
944 SDDS_CopyString(&description, s);
945
946 if (SDDS_NumberOfErrors() ||
947 (index1 = SDDS_DefineColumn(SDDSout, name, symbol, "degrees", description, NULL, SDDS_DOUBLE, 0)) < 0)
948 return 0;
949 unwrappedArgOffset = index1 - index0;
950 free(name);
951 free(symbol);
952 free(description);
953 }
954 }
955
956 free(origSymbol);
957 return 1;
958}
959
960long expandComplexColumnPairNames(SDDS_DATASET *SDDSin, char **name, char ***realName, char ***imagName, long names,
961 char **excludeName, long excludeNames, long typeMode, long typeValue) {
962 long i, j, k, realNames, imagNames, names2;
963 char **realName1, **imagName1, **realName2, **imagName2;
964 char *realPattern, *imagPattern = NULL;
965 long longest;
966
967 if (!names || !name)
968 return 0;
969 realName1 = imagName1 = realName2 = imagName2 = NULL;
970 realNames = imagNames = names2 = 0;
971 for (i = longest = 0; i < names; i++) {
972 if (strlen(name[i]) > longest)
973 longest = strlen(name[i]);
974 }
975 longest += 10;
976 if (!(realPattern = SDDS_Malloc(sizeof(*realPattern) * longest)) ||
977 !(imagPattern = SDDS_Malloc(sizeof(*imagPattern) * longest)))
978 SDDS_Bomb("Memory allocation failure");
979
980 for (i = 0; i < names; i++) {
981 for (j = 0; j < 2; j++) {
982 if (j == 0) {
983 sprintf(realPattern, "Real%s", name[i]);
984 sprintf(imagPattern, "Imag%s", name[i]);
985 } else {
986 sprintf(realPattern, "%sReal", name[i]);
987 sprintf(imagPattern, "%sImag", name[i]);
988 }
989 switch (typeMode) {
990 case FIND_ANY_TYPE:
991 case FIND_NUMERIC_TYPE:
992 case FIND_INTEGER_TYPE:
993 case FIND_FLOATING_TYPE:
994 realNames = SDDS_MatchColumns(SDDSin, &realName1, SDDS_MATCH_STRING, typeMode, realPattern, SDDS_0_PREVIOUS | SDDS_OR);
995 imagNames = SDDS_MatchColumns(SDDSin, &imagName1, SDDS_MATCH_STRING, typeMode, imagPattern, SDDS_0_PREVIOUS | SDDS_OR);
996 break;
997 case FIND_SPECIFIED_TYPE:
998 if (!SDDS_VALID_TYPE(typeValue))
999 SDDS_Bomb("Invalid type value in expandColumnPairNames");
1000 realNames = SDDS_MatchColumns(SDDSin, &realName1, SDDS_MATCH_STRING, typeMode, typeValue, realPattern, SDDS_0_PREVIOUS | SDDS_OR);
1001 imagNames = SDDS_MatchColumns(SDDSin, &imagName1, SDDS_MATCH_STRING, typeMode, typeValue, imagPattern, SDDS_0_PREVIOUS | SDDS_OR);
1002 break;
1003 default:
1004 SDDS_Bomb("Invalid typeMode in expandColumnPairNames");
1005 exit(EXIT_FAILURE);
1006 break;
1007 }
1008 if (realNames == 0)
1009 continue;
1010 if (realNames == -1 || imagNames == -1) {
1011 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
1012 SDDS_Bomb("Unable to perform column name match in expandColumnPairNames");
1013 }
1014 if (realNames != imagNames)
1015 SDDS_Bomb("Found different number of real and imaginary columns");
1016 if (excludeNames) {
1017 for (j = 0; j < excludeNames; j++)
1018 for (k = 0; k < realNames; k++)
1019 if (wild_match(realName1[k], excludeName[j])) {
1020 free(realName1[k]);
1021 free(imagName1[k]);
1022 imagName1[k] = realName1[k] = NULL;
1023 }
1024 }
1025 moveToStringArrayComplex(&realName2, &imagName2, &names2, realName1, imagName1, realNames);
1026 free(realName1);
1027 free(imagName1);
1028 }
1029 }
1030 free(realPattern);
1031 free(imagPattern);
1032 if (names2 == 0)
1033 return 0;
1034 *realName = realName2;
1035 *imagName = imagName2;
1036 return names2;
1037}
1038
1039void moveToStringArrayComplex(char ***targetReal, char ***targetImag, long *targets, char **sourceReal, char **sourceImag, long sources) {
1040 long i, j;
1041 if (!sources)
1042 return;
1043 if (!(*targetReal = SDDS_Realloc(*targetReal, sizeof(**targetReal) * (*targets + sources))) ||
1044 !(*targetImag = SDDS_Realloc(*targetImag, sizeof(**targetImag) * (*targets + sources))))
1045 SDDS_Bomb("Memory allocation failure");
1046 for (i = 0; i < sources; i++) {
1047 if (sourceReal[i] == NULL || sourceImag[i] == NULL)
1048 continue;
1049 for (j = 0; j < *targets; j++)
1050 if (strcmp(sourceReal[i], (*targetReal)[j]) == 0)
1051 break;
1052 if (j == *targets) {
1053 (*targetReal)[j] = sourceReal[i];
1054 (*targetImag)[j] = sourceImag[i];
1055 *targets += 1;
1056 }
1057 }
1058}
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_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_FreeStringArray
int32_t SDDS_FreeStringArray(char **string, int64_t strings)
Frees an array of strings by deallocating each individual string.
Definition SDDS_utils.c:2931
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
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
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
free_scanargs
void free_scanargs(SCANNED_ARG **scanned, int argc)
Definition scanargs.c:588
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