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