SDDS ToolKit Programs and Libraries for C and Python
•All Classes Files Functions Variables Macros Pages
sddsdigfilter.c
Go to the documentation of this file.
1/**
2 * @file sddsdigfilter.c
3 * @brief General series/cascade time-domain filtering utility.
4 *
5 * @details
6 * This program applies series/cascade time-domain filtering to input data using proportional,
7 * low-pass, high-pass, digital, and analog filters. It supports cascading of filters and
8 * customization for various filter configurations. Data can be supplied via files or pipes.
9 *
10 * @section Usage
11 * ```
12 * sddsdigfilter [<inputfile>] [<outputfile>]
13 * [-pipe=[input][,output]]
14 * -columns=<xcol>,<ycol>
15 * [-proportional=<gain>]
16 * [-lowpass=<gain>,<cutoff_freq>]
17 * [-highpass=<gain>,<cutoff_freq>]
18 * [-digitalfilter=<sddsfile>,<a_coeff_col>,<b_coeff_col>]
19 * [-digitalfilter=[A,<a0>,<a1>[,...]][,B,<b0>,<b1>[,...]]]
20 * [-analogfilter=<sddsfile>,<c_coeff_col>,<d_coeff_col>]
21 * [-analogfilter=[C,<c0>,<c1>[,...]][,D,<d0>,<d1>[,...]]]
22 * [-cascade]
23 * [-verbose]
24 * [-majorOrder=row|column]
25 * ```
26 *
27 * @section Options
28 * | Required | Description |
29 * |-------------------------------------|-------------------------------------------------------------|
30 * | `-columns` | Specifies the x and y columns to process. |
31 *
32 * | Optional | Description |
33 * |-------------------------------------|-------------------------------------------------------------|
34 * | `-pipe` | Uses standard input/output. |
35 * | `-proportional` | Applies a proportional gain to the input data. |
36 * | `-lowpass` | Applies a low-pass filter with specified gain and cutoff. |
37 * | `-highpass` | Applies a high-pass filter with specified gain and cutoff. |
38 * | `-digitalfilter` | Apply a digital filter with given coefficients. |
39 * | `-analogfilter` | Apply an analog filter with given coefficients. |
40 * | `-cascade` | Enables cascading of filters. |
41 * | `-verbose` | Outputs detailed processing information. |
42 * | `-majorOrder` | Specifies the major order of data (row or column). |
43 *
44 * @subsection Incompatibilities
45 * - `-cascade` is incompatible with:
46 * - Using no filters.
47 *
48 * @subsection MEO Mutually Exclusive Options
49 * - Only one of the following may be specified:
50 * - `-digitalfilter=<sddsfile>,<a_col>,<b_col>`
51 * - `-digitalfilter=[A,<a0>,<a1>,...][,B,<b0>,<b1>,...]`
52 * - Only one of the following may be specified:
53 * - `-analogfilter=<sddsfile>,<c_col>,<d_col>`
54 * - `-analogfilter=[C,<c0>,<c1>,...][,D,<d0>,<d1>,...]`
55 *
56 * @note Digital filter form:
57 * @verbatim
58 * b0 + b1.z^-1 + ... + bn.z^-n
59 * H(z) = --------------------------------
60 * a0 + a1.z^-1 + ... + an.z^-n
61 * @endverbatim
62 *
63 * @note Analog filter form:
64 * @verbatim
65 * d0 + d1.s^1 + ... + dn.s^n
66 * H(s) = ------------------------------
67 * c0 + c1.s^1 + ... + cn.s^n
68 * @endverbatim
69 *
70 * @copyright
71 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
72 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
73 *
74 * @license
75 * This file is distributed under the terms of the Software License Agreement
76 * found in the file LICENSE included with this distribution.
77 *
78 * @authors
79 * J. Carwardine, C. Saunders, R. Soliday, Xuesong Jiao, H. Shang, L. Emery
80 */
81
82#include "mdb.h"
83#include "SDDS.h"
84#include "scan.h"
85#include "SDDSutils.h"
86
87#define TRUE 1
88#define FALSE 0
89
90/* Enumeration for option types */
91enum option_type {
92 CLO_LOWPASS,
93 CLO_HIGHPASS,
94 CLO_PROPORTIONAL,
95 CLO_ANALOGFILTER,
96 CLO_DIGITALFILTER,
97 CLO_CASCADE,
98 CLO_COLUMN,
99 CLO_PIPE,
100 CLO_VERBOSE,
101 CLO_MAJOR_ORDER,
102 N_OPTIONS
103};
104
105static char *option[N_OPTIONS] = {
106 "lowpass",
107 "highpass",
108 "proportional",
109 "analogfilter",
110 "digitalfilter",
111 "cascade",
112 "columns",
113 "pipe",
114 "verbose",
115 "majorOrder"
116};
117
118char *usage =
119 "\nsddsdigfilter [<inputfile>] [<outputfile>] [-pipe=[input][,output]]\n"
120 " -columns=<xcol>,<ycol>[,...]\n"
121 " [options]:\n"
122 " -proportional=<gain>\n"
123 " -lowpass=<gain>,<cutoff_freq>\n"
124 " -highpass=<gain>,<cutoff_freq>\n"
125 " -digitalfilter=<sddsfile>,<a_coeff_col>,<b_coeff_col>\n"
126 " or -digitalfilter=[A,<a0>,<a1>,...][,B,<b0>,<b1>,...]\n"
127 " -analogfilter=<sddsfile>,<c_coeff_col>,<d_coeff_col>\n"
128 " or -analogfilter=[C,<c0>,<c1>,...][,D,<d0>,<d1>,...]\n"
129 " -cascade\n"
130 " -verbose\n"
131 " -majorOrder=row|column\n"
132 "\nDigital filter form:\n"
133 " b0 + b1.z^-1 + ... + bn.z^-n\n"
134 " --------------------------------\n"
135 " a0 + a1.z^-1 + ... + an.z^-n\n"
136 "\nAnalog filter form:\n"
137 " d0 + d1.s^1 + ... + dn.s^n\n"
138 " ------------------------------\n"
139 " c0 + c1.s^1 + ... + cn.s^n\n\n";
140
141/**************************/
142
143/* structure definitions */
144
145typedef struct
146{
147 double gain, cutoff;
148} LOWHIGH;
149
150typedef struct
151{
152 double gain;
153} GAIN;
154
155typedef struct
156{
157 double *ACcoeff;
158 double *BDcoeff;
159 int64_t ncoeffs;
160} ANADIG;
161
162typedef struct
163{
164 void **filter;
165 long *type, numFilters;
166 double *time, *input, *output;
167} STAGE;
168
169long processFilterParams(char **args, long items, char *fistStr, char *nextStr, ANADIG *ADptr);
170long processFilterParamsFromFile(char *filename, char *ACColumn, char *BDColumn, ANADIG *ADptr);
171
172/* filter functions */
173long lowpass_function(STAGE *stage, long filterNum, int64_t npoints);
174long highpass_function(STAGE *stage, long filterNum, int64_t npoints);
175long proportional_function(STAGE *stage, long filterNum, int64_t npoints);
176long analog_function(STAGE *stage, long filterNum, int64_t npoints);
177long digital_function(STAGE *stage, long filterNum, int64_t npoints);
178
179/* signal processing functions */
180long response(double *input, double *output, int64_t nInput, int64_t nOutput, double *Bcoeff, double *Acoeff, int64_t nBcoeffs, int64_t nAcoeffs);
181
182void dspbiln(double *d, double *c, int64_t ln, double *b, double *a, double *work, long *error);
183void dspfblt(double *d, double *c, long ln, long iband, double fln, double fhn, double *b, double *a, double *work, long *error);
184double dspbfct(long *i1, long *i2);
185double dpow_ri(double *ap, long *bp);
186
187long (*filter_funct[])(STAGE *stage, long filterNum, int64_t npoints) =
188 {
189 lowpass_function, highpass_function, proportional_function, analog_function, digital_function};
190
191/**************************/
192long verbose;
193
194int main(int argc, char *argv[]) {
195 SCANNED_ARG *s_arg;
196 SDDS_DATASET SDDS_input, SDDS_output;
197 char *inputfile, *outputfile, *xcolName, **ycolName, **ycolMatch;
198 long columnsupplied, totalFilterCount, totalStages, tmpfileUsed;
199 int64_t npoints;
200 int32_t yColumns, ycolMatches;
201 long i, j, i_arg, filterNum, stageNum, currentFilter, currentStage, filterType, error;
202 unsigned long pipeFlags, majorOrderFlag;
203 char outColName[256], outColDesc[256];
204 double *xcol, *ycol;
205 short columnMajorOrder = -1;
206
207 LOWHIGH *LHptr;
208 GAIN *Gptr;
209 ANADIG *ADptr;
210 STAGE *stage;
211
212 SDDS_RegisterProgramName(argv[0]);
213 argc = scanargs(&s_arg, argc, argv);
214 if (argc < 2)
215 bomb(NULL, usage);
216
217 /* initialize flags */
218 tmpfileUsed = 0;
219 pipeFlags = 0;
220 verbose = 0;
221 columnsupplied = 0;
222 inputfile = outputfile = xcolName = NULL;
223 ycolName = ycolMatch = NULL;
224 yColumns = ycolMatches = 0;
225
226 /* set up stage counters etc */
227 totalFilterCount = 0;
228 filterNum = 0;
229 stageNum = 0;
230
231 /* allocate memory for first stage */
232 stage = (STAGE *)malloc(sizeof(*stage));
233 stage[0].filter = (void *)malloc(sizeof(*stage[0].filter));
234 stage[0].type = (long *)malloc(sizeof(*stage[0].type));
235
236 for (i_arg = 1; i_arg < argc; i_arg++) {
237 if (s_arg[i_arg].arg_type == OPTION) {
238 delete_chars(s_arg[i_arg].list[0], "_");
239 switch (match_string(s_arg[i_arg].list[0], option, N_OPTIONS, 0)) {
240 case CLO_MAJOR_ORDER:
241 majorOrderFlag = 0;
242 s_arg[i_arg].n_items--;
243 if (s_arg[i_arg].n_items > 0 && (!scanItemList(&majorOrderFlag, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0, "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER, "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
244 SDDS_Bomb("invalid -majorOrder syntax/values");
245 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
246 columnMajorOrder = 1;
247 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
248 columnMajorOrder = 0;
249 break;
250 case CLO_LOWPASS:
251 if (s_arg[i_arg].n_items < 2)
252 SDDS_Bomb("Invalid -lowpass option.");
253 LHptr = (LOWHIGH *)malloc(sizeof(*LHptr));
254 if (!get_double(&(LHptr->gain), s_arg[i_arg].list[1]))
255 SDDS_Bomb("Invalid -lowpass value provided.");
256 if (s_arg[i_arg].n_items > 2 && !get_double(&(LHptr->cutoff), s_arg[i_arg].list[2]))
257 SDDS_Bomb("Invalid -lowpass value provided.");
258 stage[stageNum].filter[filterNum] = LHptr;
259 stage[stageNum].type[filterNum] = CLO_LOWPASS;
260 stage[stageNum].numFilters = ++filterNum;
261 totalFilterCount++;
262 /* allocate memory for next filter */
263 stage[stageNum].filter = (void *)realloc(stage[stageNum].filter, sizeof(*stage[stageNum].filter) * (filterNum + 1));
264 stage[stageNum].type = (long *)realloc(stage[stageNum].type, sizeof(*stage[stageNum].type) * (filterNum + 1));
265 break;
266 case CLO_HIGHPASS:
267 if (s_arg[i_arg].n_items < 2)
268 SDDS_Bomb("Invalid -highpass option.");
269 LHptr = (LOWHIGH *)malloc(sizeof(*LHptr));
270 if (!get_double(&(LHptr->gain), s_arg[i_arg].list[1]))
271 SDDS_Bomb("Invalid -highpass value provided.");
272 if (s_arg[i_arg].n_items > 2 && !get_double(&(LHptr->cutoff), s_arg[i_arg].list[2]))
273 SDDS_Bomb("Invalid -highpass value provided.");
274 stage[stageNum].filter[filterNum] = LHptr;
275 stage[stageNum].type[filterNum] = CLO_HIGHPASS;
276 stage[stageNum].numFilters = ++filterNum;
277 totalFilterCount++;
278 /* allocate memory for next filter */
279 stage[stageNum].filter = (void *)realloc(stage[stageNum].filter, sizeof(*stage[stageNum].filter) * (filterNum + 1));
280 stage[stageNum].type = (long *)realloc(stage[stageNum].type, sizeof(*stage[stageNum].type) * (filterNum + 1));
281 break;
282 case CLO_PROPORTIONAL:
283 if (s_arg[i_arg].n_items != 2)
284 SDDS_Bomb("Invalid -proportional option.");
285 Gptr = (GAIN *)malloc(sizeof(*Gptr));
286 if (!get_double(&(Gptr->gain), s_arg[i_arg].list[1]))
287 SDDS_Bomb("Invalid -proportional value provided.");
288 stage[stageNum].filter[filterNum] = Gptr;
289 stage[stageNum].type[filterNum] = CLO_PROPORTIONAL;
290 stage[stageNum].numFilters = ++filterNum;
291 totalFilterCount++;
292 /* allocate memory for next filter */
293 stage[stageNum].filter = (void *)realloc(stage[stageNum].filter, sizeof(*stage[stageNum].filter) * (filterNum + 1));
294 stage[stageNum].type = (long *)realloc(stage[stageNum].type, sizeof(*stage[stageNum].type) * (filterNum + 1));
295 break;
296 case CLO_ANALOGFILTER:
297 if (s_arg[i_arg].n_items < 2)
298 SDDS_Bomb("Invalid -analogfilter option.");
299 ADptr = (ANADIG *)malloc(sizeof(*ADptr));
300 ADptr->ACcoeff = ADptr->BDcoeff = NULL;
301 if (strcasecmp(s_arg[i_arg].list[1], "C") == 0 || strcasecmp(s_arg[i_arg].list[1], "D") == 0) {
302 processFilterParams(s_arg[i_arg].list, s_arg[i_arg].n_items, "C", "D", ADptr);
303 } else {
304 if (s_arg[i_arg].n_items != 4)
305 SDDS_Bomb("Invalid -analogfilter option for providing coefficient file and coefficient column names.");
306 processFilterParamsFromFile(s_arg[i_arg].list[1], s_arg[i_arg].list[2], s_arg[i_arg].list[3], ADptr);
307 }
308 stage[stageNum].filter[filterNum] = ADptr;
309 stage[stageNum].type[filterNum] = CLO_ANALOGFILTER;
310 stage[stageNum].numFilters = ++filterNum;
311 totalFilterCount++;
312 /* allocate memory for next filter */
313 stage[stageNum].filter = (void *)realloc(stage[stageNum].filter, sizeof(*stage[stageNum].filter) * (filterNum + 1));
314 stage[stageNum].type = (long *)realloc(stage[stageNum].type, sizeof(*stage[stageNum].type) * (filterNum + 1));
315 break;
316 case CLO_DIGITALFILTER:
317 if (s_arg[i_arg].n_items < 2)
318 SDDS_Bomb("Invalid -digitalfilter option.");
319 ADptr = (ANADIG *)malloc(sizeof(*ADptr));
320 ADptr->ACcoeff = ADptr->BDcoeff = NULL;
321 if (strcasecmp(s_arg[i_arg].list[1], "A") == 0 || strcasecmp(s_arg[i_arg].list[1], "B") == 0) {
322 processFilterParams(s_arg[i_arg].list, s_arg[i_arg].n_items, "A", "B", ADptr);
323 } else {
324 if (s_arg[i_arg].n_items != 4)
325 SDDS_Bomb("Invalid -digitalfilter option for providing coefficient file and coefficient column names.");
326 processFilterParamsFromFile(s_arg[i_arg].list[1], s_arg[i_arg].list[2], s_arg[i_arg].list[3], ADptr);
327 }
328 stage[stageNum].filter[filterNum] = ADptr;
329 stage[stageNum].type[filterNum] = CLO_DIGITALFILTER;
330 stage[stageNum].numFilters = ++filterNum;
331 totalFilterCount++;
332 /* allocate memory for next filter */
333 stage[stageNum].filter = (void *)realloc(stage[stageNum].filter, sizeof(*stage[stageNum].filter) * (filterNum + 1));
334 stage[stageNum].type = (long *)realloc(stage[stageNum].type, sizeof(*stage[stageNum].type) * (filterNum + 1));
335 break;
336 case CLO_CASCADE:
337 stageNum++;
338 /* allocate memory for next stage */
339 stage = (STAGE *)realloc(stage, sizeof(*stage) * (stageNum + 1));
340 stage[stageNum].filter = (void *)malloc(sizeof(*stage[stageNum].filter));
341 stage[stageNum].type = (long *)malloc(sizeof(*stage[stageNum].type));
342 filterNum = 0;
343 stage[stageNum].numFilters = filterNum;
344 break;
345 case CLO_COLUMN:
346 if (s_arg[i_arg].n_items < 3)
347 SDDS_Bomb("Invalid -column option.");
348 xcolName = s_arg[i_arg].list[1];
349 ycolMatches = s_arg[i_arg].n_items - 2;
350 ycolMatch = malloc(sizeof(*ycolMatch) * ycolMatches);
351 for (i = 2; i < s_arg[i_arg].n_items; i++)
352 ycolMatch[i - 2] = s_arg[i_arg].list[i];
353 columnsupplied = 1;
354 break;
355 case CLO_PIPE:
356 if (!processPipeOption(s_arg[i_arg].list + 1, s_arg[i_arg].n_items - 1, &pipeFlags))
357 SDDS_Bomb("invalid -pipe syntax");
358 break;
359 case CLO_VERBOSE:
360 verbose = 1;
361 break;
362 default:
363 fprintf(stderr, "Error (%s): unknown switch: %s\n", argv[0], s_arg[i_arg].list[0]);
364 exit(EXIT_FAILURE);
365 break;
366 }
367 } else {
368 if (inputfile == NULL)
369 inputfile = s_arg[i_arg].list[0];
370 else if (outputfile == NULL)
371 outputfile = s_arg[i_arg].list[0];
372 else
373 SDDS_Bomb("too many filenames.");
374 }
375 }
376
377 processFilenames("sddsdigfilter", &inputfile, &outputfile, pipeFlags, 1, &tmpfileUsed);
378 /* check options */
379 if (totalFilterCount == 0 || columnsupplied == 0) {
380 fprintf(stderr, "no filter or no columns supplied.\n");
381 exit(EXIT_FAILURE);
382 }
383 totalStages = stageNum + 1;
384 if (!SDDS_InitializeInput(&SDDS_input, inputfile))
385 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
386 if (SDDS_CheckColumn(&SDDS_input, xcolName, NULL, 0, stderr) != SDDS_CHECK_OKAY) {
387 fprintf(stderr, "xColumn %s does not exist.\n", xcolName);
388 exit(EXIT_FAILURE);
389 }
390 ycolName = getMatchingSDDSNames(&SDDS_input, ycolMatch, ycolMatches, &yColumns, SDDS_MATCH_COLUMN);
391 for (i = 0; i < yColumns; i++) {
392 if (SDDS_CheckColumn(&SDDS_input, ycolName[i], NULL, 0, stderr) != SDDS_CHECK_OKAY) {
393 fprintf(stderr, "yColumn %s does not exist.\n", ycolName[i]);
394 exit(EXIT_FAILURE);
395 }
396 }
397 if (!SDDS_InitializeCopy(&SDDS_output, &SDDS_input, outputfile, "w"))
398 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
399 if (columnMajorOrder != -1)
400 SDDS_output.layout.data_mode.column_major = columnMajorOrder;
401 else
402 SDDS_output.layout.data_mode.column_major = SDDS_input.layout.data_mode.column_major;
403 for (i = 0; i < yColumns; i++) {
404 sprintf(outColName, "DigFiltered%s", ycolName[i]);
405 sprintf(outColDesc, "Digital Filtered %s", ycolName[i]);
406 if (SDDS_DefineColumn(&SDDS_output, outColName, NULL, NULL, outColDesc, NULL, SDDS_DOUBLE, 0) < 0)
407 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
408 }
409 if (!SDDS_WriteLayout(&SDDS_output))
410 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
411
412 while ((SDDS_ReadPage(&SDDS_input)) > 0) {
413 if (!SDDS_CopyPage(&SDDS_output, &SDDS_input))
414 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
415 if ((npoints = SDDS_CountRowsOfInterest(&SDDS_input)) < 0)
416 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
417 if (npoints) {
418 /**** get input data ****/
419 if (!(xcol = SDDS_GetColumnInDoubles(&SDDS_input, xcolName)))
420 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
421 for (i = 0; i < yColumns; i++) {
422 sprintf(outColName, "DigFiltered%s", ycolName[i]);
423 if (!(ycol = SDDS_GetColumnInDoubles(&SDDS_input, ycolName[i])))
424 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
425 /**** PROCESS EACH STAGE ****/
426 for (currentStage = 0; currentStage < totalStages; currentStage++) {
427 if (verbose == 1) {
428 fprintf(stdout, "STAGE %ld\n", currentStage);
429 }
430 /* set up time array */
431 stage[currentStage].time = xcol;
432 /* set up input data for current stage */
433 if (currentStage == 0) {
434 stage[currentStage].input = ycol;
435 } else {
436 stage[currentStage].input = stage[currentStage - 1].output;
437 }
438 /* allocate memory for output array */
439 stage[currentStage].output = (double *)calloc(npoints, sizeof(*stage[currentStage].output));
440 /**** PROCESS EACH FILTER IN CURRENT STAGE ****/
441 for (currentFilter = 0; currentFilter < stage[currentStage].numFilters; currentFilter++) {
442 if (verbose == 1) {
443 fprintf(stdout, "filter %ld\n", currentFilter);
444 }
445 /* send data to relevant filter function */
446 filterType = stage[currentStage].type[currentFilter];
447 error = (*filter_funct[filterType])(&(stage[currentStage]), currentFilter, npoints);
448 if (error != 0) {
449 fprintf(stderr, "nyquist violation.\n");
450 exit(EXIT_FAILURE);
451 }
452 /* end of filter loop */
453 }
454 /* end of stage loop */
455 }
456 /**** End of Processing ****/
457 /* produce outputfile */
458 if (!SDDS_SetColumnFromDoubles(&SDDS_output, SDDS_BY_NAME, stage[totalStages - 1].output, npoints, outColName))
459 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
460 /* clean up */
461 free(ycol);
462 for (currentStage = 0; currentStage < totalStages; currentStage++)
463 free(stage[currentStage].output);
464 }
465 free(xcol);
466 }
467 if (!SDDS_WritePage(&SDDS_output))
468 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
469 }
470 if (!SDDS_Terminate(&SDDS_input) || !SDDS_Terminate(&SDDS_output))
471 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
472 for (i = 0; i < totalStages; i++) {
473 for (j = 0; j < stage[i].numFilters; j++) {
474 switch (stage[i].type[j]) {
475 case CLO_ANALOGFILTER:
476 case CLO_DIGITALFILTER:
477 free(((ANADIG *)(stage[i].filter[j]))->ACcoeff);
478 free(((ANADIG *)(stage[i].filter[j]))->BDcoeff);
479 free((ANADIG *)(stage[i].filter[j]));
480 break;
481 case CLO_PROPORTIONAL:
482 free((GAIN *)(stage[i].filter[j]));
483 break;
484 case CLO_HIGHPASS:
485 case CLO_LOWPASS:
486 free((LOWHIGH *)(stage[i].filter[j]));
487 break;
488 }
489 }
490 free(stage[i].filter);
491 free(stage[i].type);
492 }
493 free(stage);
494 SDDS_FreeStringArray(ycolName, yColumns);
495 free(ycolName);
496 free(ycolMatch);
497 free_scanargs(&s_arg, argc);
498 return (EXIT_SUCCESS);
499 /**** end of main ****/
500}
501
502/********************************************/
503long lowpass_function(STAGE *stage, long filterNum, int64_t npoints) {
504
505 double C_coeff[2] = {1.0, 1.0};
506 double D_coeff[2] = {1.0, 0.0};
507
508 double A_coeff[2], B_coeff[2];
509 double timePerPoint, normFrequency;
510 double *work, *tmpoutput;
511
512 LOWHIGH *filter_ptr;
513
514 long error;
515 int64_t i;
516
517 filter_ptr = stage->filter[filterNum];
518
519 timePerPoint = (stage->time[1]) - (stage->time[0]);
520 if (timePerPoint < 0)
521 timePerPoint *= -1.0;
522
523 if (filter_ptr->gain == 0.0) {
524 for (i = 0; i < npoints; i++) {
525 stage->output[i] = 0.0;
526 }
527 return 0;
528 }
529
530 work = (double *)calloc(npoints, sizeof(*work));
531 tmpoutput = (double *)calloc(npoints, sizeof(*tmpoutput));
532
533 normFrequency = timePerPoint * (filter_ptr->cutoff);
534
535 /* create digital filter coeffs */
536 dspfblt(D_coeff, C_coeff, 1, 1, normFrequency, 0.0, B_coeff, A_coeff, work, &error);
537 if (error != 0) {
538 free(work);
539 free(tmpoutput);
540 return error;
541 }
542 /* shift A coeffs by one */
543 A_coeff[1] = A_coeff[0];
544 A_coeff[0] = 1.0;
545
546 if (verbose == 1) {
547 fprintf(stdout, "a0: %f\tb0: %f\n", A_coeff[0], B_coeff[0]);
548 fprintf(stdout, "a1: %f\tb1: %f\n\n", A_coeff[1], B_coeff[1]);
549 }
550
551 /* apply digital filter to data */
552 if (response(stage->input, tmpoutput, npoints, npoints, B_coeff, A_coeff, 2, 2)) {
553 fprintf(stderr, "Invalid coeffs in lowpass filter.\n");
554 free(work);
555 free(tmpoutput);
556 exit(EXIT_FAILURE);
557 }
558
559 /* sum this output with the existing stage output */
560 for (i = 0; i < npoints; i++) {
561 stage->output[i] += (filter_ptr->gain) * tmpoutput[i];
562 }
563
564 /* clean up */
565 free(work);
566 free(tmpoutput);
567
568 return 0;
569 /* end of lowpass_function() */
570}
571
572/************************************************/
573long highpass_function(STAGE *stage, long filterNum, int64_t npoints) {
574
575 double C_coeff[] = {1.0, 1.0};
576 double D_coeff[] = {1.0, 0.0};
577
578 double A_coeff[2], B_coeff[2];
579 double timePerPoint, normFrequency;
580 double *work, *tmpoutput;
581
582 LOWHIGH *filter_ptr;
583
584 long error;
585 int64_t i;
586
587 filter_ptr = stage->filter[filterNum];
588
589 timePerPoint = (stage->time[1]) - (stage->time[0]);
590 if (timePerPoint < 0)
591 timePerPoint *= -1.0;
592
593 if (filter_ptr->gain == 0.0) {
594 for (i = 0; i < npoints; i++) {
595 stage->output[i] = 0.0;
596 }
597 return 0;
598 }
599
600 work = (double *)calloc(npoints, sizeof(*work));
601 tmpoutput = (double *)calloc(npoints, sizeof(*tmpoutput));
602
603 normFrequency = timePerPoint * (filter_ptr->cutoff);
604
605 /* create digital filter coeffs */
606 dspfblt(D_coeff, C_coeff, 1, 2, normFrequency, 0.0, B_coeff, A_coeff, work, &error);
607 if (error != 0) {
608 free(work);
609 free(tmpoutput);
610 return error;
611 }
612 /* shift A coeffs by one */
613 A_coeff[1] = A_coeff[0];
614 A_coeff[0] = 1.0;
615
616 if (verbose == 1) {
617 fprintf(stdout, "a0: %6.6g\tb0: %6.6g\n", A_coeff[0], B_coeff[0]);
618 fprintf(stdout, "a1: %6.6g\tb1: %6.6g\n\n", A_coeff[1], B_coeff[1]);
619 }
620
621 /* apply digital filter to data */
622 if (response(stage->input, tmpoutput, npoints, npoints, B_coeff, A_coeff, 2, 2)) {
623 fprintf(stderr, "invalid coeffs in highpass filter.\n");
624 free(work);
625 free(tmpoutput);
626 exit(EXIT_FAILURE);
627 }
628
629 /* sum this output with the existing stage output */
630 for (i = 0; i < npoints; i++) {
631 stage->output[i] += (filter_ptr->gain) * tmpoutput[i];
632 }
633
634 /* clean up */
635 free(work);
636 free(tmpoutput);
637
638 return 0;
639 /* end of highpass_function() */
640}
641
642/*****************************************/
643
644long digital_function(STAGE *stage, long filterNum, int64_t npoints) {
645
646 double *tmpoutput;
647
648 ANADIG *filter_ptr;
649
650 int64_t i;
651
652 filter_ptr = stage->filter[filterNum];
653
654 tmpoutput = (double *)calloc(npoints, sizeof(*tmpoutput));
655
656 if (verbose == 1) {
657 for (i = 0; i < filter_ptr->ncoeffs; i++) {
658 fprintf(stdout, "a%" PRId64 ": %6.6g\tb%" PRId64 " %6.6g\n", i, filter_ptr->ACcoeff[i], i, filter_ptr->BDcoeff[i]);
659 }
660 }
661
662 /* apply digital filter to data */
663 if (response(stage->input, tmpoutput, npoints, npoints, filter_ptr->BDcoeff, filter_ptr->ACcoeff, (filter_ptr->ncoeffs), (filter_ptr->ncoeffs))) {
664 fprintf(stderr, "Invalid coefficients in digitalfilter.\n");
665 free(tmpoutput);
666 exit(EXIT_FAILURE);
667 }
668
669 /* sum this output with the existing stage output */
670 for (i = 0; i < npoints; i++) {
671 stage->output[i] += tmpoutput[i];
672 }
673
674 /* clean up */
675 free(tmpoutput);
676
677 return 0;
678 /* end of digital_function() */
679}
680
681/*****************************************/
682
683long analog_function(STAGE *stage, long filterNum, int64_t npoints) {
684 double sampleFrequency, sampleTime;
685 double *tmpoutput, *work;
686 double *Bcoeff, *Acoeff;
687 ANADIG *filter_ptr;
688
689 long error;
690 int64_t i;
691
692 filter_ptr = stage->filter[filterNum];
693 work = (double *)calloc(npoints, sizeof(*work));
694 tmpoutput = (double *)calloc(npoints, sizeof(*tmpoutput));
695 Bcoeff = (double *)calloc(filter_ptr->ncoeffs, sizeof(*Bcoeff));
696 Acoeff = (double *)calloc(filter_ptr->ncoeffs, sizeof(*Acoeff));
697
698 /* normalize analog frequency components */
699 /* this is needed because the bilinear transform routine
700 does not take into account the T/2 in the transformation */
701 sampleTime = (stage->time[1] - stage->time[0]);
702 sampleFrequency = 2.0 / sampleTime;
703 if (sampleFrequency < 0.0)
704 sampleFrequency *= -1.0;
705
706 for (i = 0; i < filter_ptr->ncoeffs; i++) {
707 filter_ptr->BDcoeff[i] *= pow(sampleFrequency, i);
708 filter_ptr->ACcoeff[i] *= pow(sampleFrequency, i);
709 }
710
711 /* do bilinear transform */
712 dspbiln(filter_ptr->BDcoeff, filter_ptr->ACcoeff, (filter_ptr->ncoeffs) - 1, Bcoeff, Acoeff, work, &error);
713 if (error) {
714 fprintf(stderr, "Error(%ld): invalid or all-zero analog transfer function\n", error);
715 free(work);
716 free(tmpoutput);
717 free(Bcoeff);
718 free(Acoeff);
719 exit(EXIT_FAILURE);
720 }
721
722 /* shift A coeffs back by one (dspbiln assumes A0 is actually A1) */
723 for (i = filter_ptr->ncoeffs - 1; i > 0; i--) {
724 Acoeff[i] = Acoeff[i - 1];
725 }
726 Acoeff[0] = 1.0;
727
728 if (verbose == 1) {
729 for (i = 0; i < filter_ptr->ncoeffs; i++) {
730 fprintf(stdout, "c%" PRId64 ": %6.6g\td%" PRId64 ": %6.6g\ta%" PRId64 ": %6.6g\tb%" PRId64 ": %6.6g\n",
731 i, filter_ptr->ACcoeff[i],
732 i, filter_ptr->BDcoeff[i],
733 i, Acoeff[i],
734 i, Bcoeff[i]);
735 }
736 }
737
738 /* apply digital filter to data */
739 if (response(stage->input, tmpoutput, npoints, npoints, Bcoeff, Acoeff, (filter_ptr->ncoeffs), (filter_ptr->ncoeffs))) {
740 fprintf(stderr, "Invalid coefficients in analogfilter.\n");
741 free(work);
742 free(tmpoutput);
743 free(Bcoeff);
744 free(Acoeff);
745 exit(EXIT_FAILURE);
746 }
747
748 /* sum this output with the existing stage output */
749 for (i = 0; i < npoints; i++) {
750 stage->output[i] += tmpoutput[i];
751 }
752
753 /* clean up */
754 free(tmpoutput);
755 free(Bcoeff);
756 free(Acoeff);
757 free(work);
758
759 return 0;
760
761 /* end of analog_function() */
762}
763
764/*****************************************/
765
766long proportional_function(STAGE *stage, long filterNum, int64_t npoints) {
767
768 GAIN *filter_ptr;
769 int64_t i;
770 double *tmpoutput;
771
772 tmpoutput = (double *)calloc(npoints, sizeof(*tmpoutput));
773 filter_ptr = stage->filter[filterNum];
774
775 for (i = 0; i < npoints; i++) {
776 tmpoutput[i] = (stage->input[i]) * (filter_ptr->gain);
777 }
778
779 /* sum this output with the existing stage output */
780 for (i = 0; i < npoints; i++) {
781 stage->output[i] += tmpoutput[i];
782 }
783
784 free(tmpoutput);
785
786 return 0;
787 /* end of proportional_function() */
788}
789
790/**************************************/
791
792long response(double *input, double *output, int64_t nInput, int64_t nOutput, double *Bcoeff, double *Acoeff, int64_t nBcoeffs, int64_t nAcoeffs) {
793
794 /* This routine adapted from Stearns & David, "Digital Signal Processing
795 Algorithms in Fortran & C, Prentice Hall, 1993 */
796
797 /* realizes the equation: */
798 /* Acoeff[0] * output[k] = SUM( Bcoeff[n] * input[k] ) - SUM( Acoeff[m] * output[m] ) */
799 /* n m */
800 /* where 0 < n < nBcoeffs, and 1 < m < nAcoeffs */
801 /* normal situation is where Acoeff[0] = 1.0 */
802 /* if there are fewer points in the input than the output, */
803 /* the last point of input is continued on */
804
805 int64_t k, n, m, max, inputPoint;
806
807 if (Acoeff[0] == 0.0 || nAcoeffs < 1 || nBcoeffs < 1) {
808 return 1;
809 }
810
811 for (k = 0; k < nOutput; k++) {
812 output[k] = 0.0;
813
814 max = (k < nBcoeffs) ? k : nBcoeffs;
815 for (n = 0; n < max; n++) {
816 if ((k - n) >= 0) {
817 /* take the last supplied input point or earlier */
818 inputPoint = ((k - n) < nInput) ? (k - n) : (nInput - 1);
819 output[k] += Bcoeff[n] * input[inputPoint];
820 }
821 }
822
823 max = (k < nAcoeffs) ? k : nAcoeffs;
824 for (m = 1; m < max; m++) {
825 if ((k - m) >= 0) {
826 output[k] -= Acoeff[m] * output[k - m];
827 }
828 }
829
830 output[k] /= Acoeff[0];
831 }
832
833 return 0;
834
835 /* end of response */
836}
837
838/****************************************/
839
840/* DSPBILN 11/13/85 */
841/* This routine from Stearns & David, "Digital Signal Processing Algorithms
842 in Fortran & C, Prentice Hall, 1993 */
843/* CONVERTS ANALOG H(S) TO DIGITAL H(Z) VIA BILINEAR TRANSFORM */
844/* ANALOG TRANSFER FUNCTION DIGITAL TRANSFER FUNCTION */
845/* D(L)*S**L+.....+D(0) B(0)+......+B(L)*Z**-L */
846/* H(S)=--------------------- H(Z)=---------------------- */
847/* C(L)*S**L+.....+C(0) 1+.......+A(L)*Z**-L */
848/* H(S) IS ASSUMED TO BE PRE-SCALED AND PRE-WARPED */
849/* LN DSPECIFIES THE LENGTH OF THE COEFFICIENT ARRAYS */
850/* FILTER ORDER L IS COMPUTED INTERNALLY */
851/* WORK IS AN INTERNAL ARRAY (2D) SIZED TO MATCH COEF ARRAYS */
852/* IERROR=0 NO ERRORS DETECTED IN TRANSFORMATION */
853/* 1 ALL ZERO TRANSFER FUNCTION */
854/* 2 INVALID TRANSFER FUNCTION; Y(K) COEF=0 */
855
856void dspbiln(double *d, double *c, int64_t ln, double *b, double *a, double *work, long *error) {
857 /* local variables */
858 long zero_func;
859 int64_t i, j, l, work_dim1;
860 double scale, tmp;
861 double atmp;
862
863 scale = 0;
864
865 zero_func = TRUE;
866 for (i = ln; i >= 0 && zero_func; --i) {
867 if (c[i] != 0.0 || d[i] != 0.0) {
868 zero_func = FALSE;
869 }
870 }
871
872 if (zero_func) {
873 *error = 1;
874 return;
875 }
876
877 work_dim1 = ln + 1;
878
879 l = i + 1;
880 for (j = 0; j <= l; ++j) {
881 work[j * work_dim1] = 1.0;
882 }
883
884 tmp = 1.0;
885 for (i = 1; i <= l; ++i) {
886 tmp = tmp * (double)(l - i + 1) / (double)i;
887 work[i] = tmp;
888 }
889
890 for (i = 1; i <= l; ++i) {
891 for (j = 1; j <= l; ++j) {
892 work[i + j * work_dim1] = work[i + (j - 1) * work_dim1] - work[i - 1 + j * work_dim1] - work[i - 1 + (j - 1) * work_dim1];
893 }
894 }
895
896 for (i = l; i >= 0; --i) {
897 b[i] = 0.0;
898 atmp = 0.0;
899
900 for (j = 0; j <= l; ++j) {
901 b[i] += work[i + j * work_dim1] * d[j];
902 atmp += (double)work[i + j * work_dim1] * c[j];
903 }
904
905 scale = (double)atmp;
906
907 if (i != 0) {
908 a[i - 1] = (double)atmp;
909 }
910 }
911
912 if (scale == 0.0) {
913 *error = 2;
914 return;
915 }
916
917 b[0] /= scale;
918 for (i = 1; i <= l; ++i) {
919 b[i] /= scale;
920 a[i - 1] /= scale;
921 }
922
923 if (l < ln) {
924 for (i = l + 1; i <= ln; ++i) {
925 b[i] = 0.0;
926 a[i - 1] = 0.0;
927 }
928 }
929
930 *error = 0;
931 return;
932} /* dspbiln */
933
934/*********************************************/
935
936/* DSPFBLT 10/26/90 */
937/* This routine from Stearns & David, "Digital Signal Processing Algorithms
938 in Fortran & C, Prentice Hall, 1993 */
939/* CONVERTS NORMALIZED LP ANALOG H(S) TO DIGITAL H(Z) */
940/* ANALOG TRANSFER FUNCTION DIGITAL TRANSFER FUNCTION */
941/* D(M)*S**M+.....+D(0) B(0)+.....+B(L)*Z**-L */
942/* H(S)=-------------------- H(Z)=-------------------- */
943/* C(M)*S**M+.....+C(0) 1+......+A(L)*Z**-L */
944/* FILTER ORDER L IS COMPUTED INTERNALLY */
945/* IBAND=1 LOWPASS FLN=NORMALIZED CUTOFF IN HZ-SEC */
946/* 2 HIGHPASS FLN=NORMALIZED CUTOFF IN HZ-SEC */
947/* 3 BANDPASS FLN=LOW CUTOFF; FHN=HIGH CUTOFF */
948/* 4 BANDSTOP FLN=LOW CUTOFF; FHN=HIGH CUTOFF */
949/* LN DSPECIFIES COEFFICIENT ARRAY SIZE */
950/* WORK(0:LN,0:LN) IS A WORK ARRAY USED INTERNALLY */
951/* RETURN IERROR=0 NO ERRORS DETECTED */
952/* 1 ALL ZERO TRANSFER FUNCTION */
953/* 2 DSPBILN: INVALID TRANSFER FUNCTION */
954/* 3 FILTER ORDER EXCEEDS ARRAY SIZE */
955/* 4 INVALID FILTER TYPE PARAMETER (IBAND) */
956/* 5 INVALID CUTOFF FREQUENCY */
957
958void dspfblt(double *d, double *c, long ln, long iband, double fln, double fhn, double *b, double *a, double *work, long *error) {
959
960 long work_dim1, tmp_int, i, k, l, m, zero_func, ll, mm, ls;
961 double tmpc, tmpd, w = 0, w1, w2, w02 = 0, tmp;
962
963 if (iband < 1 || iband > 4) {
964 *error = 4;
965 return;
966 }
967 if ((fln <= 0.0 || fln > 0.5) || (iband >= 3 && (fln >= fhn || fhn > 0.5))) {
968 *error = 5;
969 return;
970 }
971
972 zero_func = TRUE;
973 for (i = ln; i >= 0 && zero_func; --i) {
974 if (c[i] != 0.0 || d[i] != 0.0) {
975 zero_func = FALSE;
976 }
977 }
978
979 if (zero_func) {
980 *error = 1;
981 return;
982 }
983
984 work_dim1 = ln + 1;
985
986 m = i + 1;
987 w1 = (double)tan(PI * fln);
988 l = m;
989 if (iband > 2) {
990 l = m * 2;
991 w2 = (double)tan(PI * fhn);
992 w = w2 - w1;
993 w02 = w1 * w2;
994 }
995
996 if (l > ln) {
997 *error = 3;
998 return;
999 }
1000
1001 switch (iband) {
1002 case 1:
1003
1004 /* SCALING S/W1 FOR LOWPASS,HIGHPASS */
1005
1006 for (mm = 0; mm <= m; ++mm) {
1007 d[mm] /= dpow_ri(&w1, &mm);
1008 c[mm] /= dpow_ri(&w1, &mm);
1009 }
1010
1011 break;
1012
1013 case 2:
1014
1015 /* SUBSTITUTION OF 1/S TO GENERATE HIGHPASS (HP,BS) */
1016
1017 for (mm = 0; mm <= (m / 2); ++mm) {
1018 tmp = d[mm];
1019 d[mm] = d[m - mm];
1020 d[m - mm] = tmp;
1021 tmp = c[mm];
1022 c[mm] = c[m - mm];
1023 c[m - mm] = tmp;
1024 }
1025
1026 /* SCALING S/W1 FOR LOWPASS,HIGHPASS */
1027
1028 for (mm = 0; mm <= m; ++mm) {
1029 d[mm] /= dpow_ri(&w1, &mm);
1030 c[mm] /= dpow_ri(&w1, &mm);
1031 }
1032
1033 break;
1034
1035 case 3:
1036
1037 /* SUBSTITUTION OF (S**2+W0**2)/(W*S) BANDPASS,BANDSTOP */
1038
1039 for (ll = 0; ll <= l; ++ll) {
1040 work[ll] = 0.0;
1041 work[ll + work_dim1] = 0.0;
1042 }
1043
1044 for (mm = 0; mm <= m; ++mm) {
1045 tmp_int = m - mm;
1046 tmpd = d[mm] * dpow_ri(&w, &tmp_int);
1047 tmpc = c[mm] * dpow_ri(&w, &tmp_int);
1048
1049 for (k = 0; k <= mm; ++k) {
1050 ls = m + mm - (k * 2);
1051 tmp_int = mm - k;
1052 tmp = dspbfct(&mm, &mm) / (dspbfct(&k, &k) * dspbfct(&tmp_int, &tmp_int));
1053 work[ls] += tmpd * dpow_ri(&w02, &k) * tmp;
1054 work[ls + work_dim1] += tmpc * dpow_ri(&w02, &k) * tmp;
1055 }
1056 }
1057
1058 for (ll = 0; ll <= l; ++ll) {
1059 d[ll] = work[ll];
1060 c[ll] = work[ll + work_dim1];
1061 }
1062
1063 break;
1064
1065 case 4:
1066
1067 /* SUBSTITUTION OF 1/S TO GENERATE HIGHPASS (HP,BS) */
1068
1069 for (mm = 0; mm <= (m / 2); ++mm) {
1070 tmp = d[mm];
1071 d[mm] = d[m - mm];
1072 d[m - mm] = tmp;
1073 tmp = c[mm];
1074 c[mm] = c[m - mm];
1075 c[m - mm] = tmp;
1076 }
1077
1078 /* SUBSTITUTION OF (S**2+W0**2)/(W*S) BANDPASS,BANDSTOP */
1079
1080 for (ll = 0; ll <= l; ++ll) {
1081 work[ll] = 0.0;
1082 work[ll + work_dim1] = 0.0;
1083 }
1084
1085 for (mm = 0; mm <= m; ++mm) {
1086 tmp_int = m - mm;
1087 tmpd = d[mm] * dpow_ri(&w, &tmp_int);
1088 tmpc = c[mm] * dpow_ri(&w, &tmp_int);
1089
1090 for (k = 0; k <= mm; ++k) {
1091 ls = m + mm - (k * 2);
1092 tmp_int = mm - k;
1093 tmp = dspbfct(&mm, &mm) / (dspbfct(&k, &k) * dspbfct(&tmp_int, &tmp_int));
1094 work[ls] += tmpd * dpow_ri(&w02, &k) * tmp;
1095 work[ls + work_dim1] += tmpc * dpow_ri(&w02, &k) * tmp;
1096 }
1097 }
1098
1099 for (ll = 0; ll <= l; ++ll) {
1100 d[ll] = work[ll];
1101 c[ll] = work[ll + work_dim1];
1102 }
1103
1104 break;
1105 }
1106
1107 dspbiln(d, c, ln, b, a, work, error);
1108
1109 return;
1110} /* dspfblt */
1111
1112/******************************************/
1113
1114/* DSPBFCT 11/14/85 */
1115/* This routine from Stearns & David, "Digital Signal Processing Algorithms
1116 in Fortran & C, Prentice Hall, 1993 */
1117/* GENERATES (I1)!/(I1-I2)!=I1*(I1-1)*...*(I1-I2+1). */
1118/* NOTE: 0!=1 AND DSPBFCT(I,I)=DSPBFCT(I,I-1)=I!. */
1119
1120double dspbfct(long *i1, long *i2) {
1121 /* Local variables */
1122 long i;
1123 double ret_val;
1124
1125 ret_val = 0.0;
1126 if (*i1 < 0 || *i2 < 0 || *i2 > *i1) {
1127 return (ret_val);
1128 }
1129
1130 ret_val = 1.0;
1131 for (i = *i1; i >= (*i1 - *i2 + 1); --i) {
1132 ret_val *= (double)i;
1133 }
1134
1135 return (ret_val);
1136
1137} /* DSPBFCT */
1138
1139/********************************************/
1140
1141double dpow_ri(double *ap, long *bp) {
1142 /* This routine from Stearns & David, "Digital Signal Processing Algorithms
1143 in Fortran & C, Prentice Hall, 1993 */
1144
1145 double pow, x;
1146 long n;
1147
1148 pow = 1;
1149 x = *ap;
1150 n = *bp;
1151
1152 if (n != 0) {
1153 if (n < 0) {
1154 if (x == 0) {
1155 return (pow);
1156 }
1157 n = -n;
1158 x = 1 / x;
1159 }
1160
1161 for (;;) {
1162 if (n & 01)
1163 pow *= x;
1164
1165 if (n >>= 1)
1166 x *= x;
1167 else
1168 break;
1169 }
1170 }
1171
1172 return (pow);
1173}
1174
1175/********************************************/
1176
1177long processFilterParams(char **args, long items, char *firstStr, char *nextStr, ANADIG *ADptr) {
1178 long a, b, aFlag, bFlag, i;
1179
1180 a = b = aFlag = bFlag = 0;
1181 ADptr->ACcoeff = ADptr->BDcoeff = NULL;
1182 for (i = 1; i < items; i++) {
1183 if (strcasecmp(args[i], firstStr) == 0) {
1184 aFlag = 1;
1185 bFlag = 0;
1186 } else if (strcasecmp(args[i], nextStr) == 0) {
1187 bFlag = 1;
1188 aFlag = 0;
1189 } else {
1190 if (aFlag) {
1191 ADptr->ACcoeff = SDDS_Realloc(ADptr->ACcoeff, sizeof(*(ADptr->ACcoeff)) * (a + 1));
1192 if (!get_double(&(ADptr->ACcoeff[a]), args[i]))
1193 SDDS_Bomb("Invalid filter coefficient provided.");
1194 a++;
1195 }
1196 if (bFlag) {
1197 ADptr->BDcoeff = SDDS_Realloc(ADptr->BDcoeff, sizeof(*(ADptr->BDcoeff)) * (b + 1));
1198 if (!get_double(&(ADptr->BDcoeff[b]), args[i]))
1199 SDDS_Bomb("Invalid filter coefficient provided.");
1200 b++;
1201 }
1202 }
1203 }
1204 if (a == 0 && b == 0)
1205 SDDS_Bomb("No coefficients provided.");
1206 if (a == 0) {
1207 ADptr->ACcoeff = SDDS_Realloc(ADptr->ACcoeff, sizeof(*(ADptr->ACcoeff)) * (a + 1));
1208 ADptr->ACcoeff[0] = 1.0;
1209 a = 1;
1210 }
1211 if (b == 0) {
1212 ADptr->BDcoeff = SDDS_Realloc(ADptr->BDcoeff, sizeof(*(ADptr->BDcoeff)) * (b + 1));
1213 ADptr->BDcoeff[0] = 1.0;
1214 b = 1;
1215 }
1216 ADptr->ncoeffs = a > b ? a : b;
1217 if (a < ADptr->ncoeffs) {
1218 ADptr->ACcoeff = SDDS_Realloc(ADptr->ACcoeff, sizeof(*(ADptr->ACcoeff)) * ADptr->ncoeffs);
1219 for (i = a; i < ADptr->ncoeffs; i++)
1220 ADptr->ACcoeff[i] = 0.0;
1221 }
1222 if (b < ADptr->ncoeffs) {
1223 ADptr->BDcoeff = SDDS_Realloc(ADptr->BDcoeff, sizeof(*(ADptr->BDcoeff)) * ADptr->ncoeffs);
1224 for (i = b; i < ADptr->ncoeffs; i++)
1225 ADptr->BDcoeff[i] = 0.0;
1226 }
1227 return 1;
1228}
1229
1230long processFilterParamsFromFile(char *filename, char *ACColumn, char *BDColumn, ANADIG *ADptr) {
1231 int64_t rows = 0;
1232 SDDS_DATASET SDDSdata;
1233 if (!SDDS_InitializeInput(&SDDSdata, filename)) {
1234 fprintf(stderr, "error: problem reading %s\n", filename);
1235 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1236 }
1237 if (SDDS_ReadPage(&SDDSdata) < 0)
1238 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1239 if (!(rows = SDDS_CountRowsOfInterest(&SDDSdata))) {
1240 fprintf(stderr, "error: problem counting rows in file %s\n", filename);
1241 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1242 }
1243 ADptr->ncoeffs = rows;
1244 ADptr->ACcoeff = ADptr->BDcoeff = NULL;
1245 if (!(ADptr->ACcoeff = SDDS_GetColumnInDoubles(&SDDSdata, ACColumn))) {
1246 fprintf(stderr, "error: unable to read AC coefficient column.\n");
1247 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1248 }
1249 if (!(ADptr->BDcoeff = SDDS_GetColumnInDoubles(&SDDSdata, BDColumn))) {
1250 fprintf(stderr, "error: unable to read BD coefficient column.\n");
1251 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1252 }
1253 if (!SDDS_Terminate(&SDDSdata))
1254 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1255 return 1;
1256}
1257
1258/************************************/
SDDS.h
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
SDDS_InitializeCopy
int32_t SDDS_InitializeCopy(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source, char *filename, char *filemode)
Definition SDDS_copy.c:40
SDDS_CopyPage
int32_t SDDS_CopyPage(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
Definition SDDS_copy.c:578
SDDS_SetColumnFromDoubles
int32_t SDDS_SetColumnFromDoubles(SDDS_DATASET *SDDS_dataset, int32_t mode, double *data, int64_t rows,...)
Sets the values for a single data column using double-precision floating-point numbers.
Definition SDDS_dataprep.c:1642
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_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:49
SDDS_Terminate
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1448
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1006
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1228
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:1715
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:899
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:2865
getMatchingSDDSNames
char ** getMatchingSDDSNames(SDDS_DATASET *dataset, char **matchName, int32_t matches, int32_t *names, short type)
Retrieves an array of matching SDDS entity names based on specified criteria.
Definition SDDS_utils.c:5556
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:4562
SDDS_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:432
SDDS_RegisterProgramName
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:288
SDDS_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:342
SDDS_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:677
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.
bomb
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
Definition bomb.c:26
get_double
int get_double(double *dptr, char *s)
Parses a double value from the given string.
Definition data_scan.c:40
delete_chars
char * delete_chars(char *s, char *t)
Removes all occurrences of characters found in string t from string s.
Definition delete_chars.c:28
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:356
processFilenames
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:390
free_scanargs
void free_scanargs(SCANNED_ARG **scanned, int argc)
Definition scanargs.c:584
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
Generated by doxygen 1.12.0