SDDSlib
Loading...
Searching...
No Matches
Classes | Macros | Enumerations | Functions | Variables
sddsdigfilter.c File Reference

General series/cascade time-domain filtering utility. More...

#include "mdb.h"
#include "SDDS.h"
#include "scan.h"
#include "SDDSutils.h"

Go to the source code of this file.

Classes

struct  LOWHIGH
 
struct  GAIN
 
struct  ANADIG
 
struct  STAGE
 

Macros

#define TRUE   1
 
#define FALSE   0
 

Enumerations

enum  option_type {
  CLO_LOWPASS , CLO_HIGHPASS , CLO_PROPORTIONAL , CLO_ANALOGFILTER ,
  CLO_DIGITALFILTER , CLO_CASCADE , CLO_COLUMN , CLO_PIPE ,
  CLO_VERBOSE , CLO_MAJOR_ORDER , N_OPTIONS
}
 

Functions

long processFilterParams (char **args, long items, char *fistStr, char *nextStr, ANADIG *ADptr)
 
long processFilterParamsFromFile (char *filename, char *ACColumn, char *BDColumn, ANADIG *ADptr)
 
long lowpass_function (STAGE *stage, long filterNum, int64_t npoints)
 
long highpass_function (STAGE *stage, long filterNum, int64_t npoints)
 
long proportional_function (STAGE *stage, long filterNum, int64_t npoints)
 
long analog_function (STAGE *stage, long filterNum, int64_t npoints)
 
long digital_function (STAGE *stage, long filterNum, int64_t npoints)
 
long response (double *input, double *output, int64_t nInput, int64_t nOutput, double *Bcoeff, double *Acoeff, int64_t nBcoeffs, int64_t nAcoeffs)
 
void dspbiln (double *d, double *c, int64_t ln, double *b, double *a, double *work, long *error)
 
void dspfblt (double *d, double *c, long ln, long iband, double fln, double fhn, double *b, double *a, double *work, long *error)
 
double dspbfct (long *i1, long *i2)
 
double dpow_ri (double *ap, long *bp)
 
int main (int argc, char *argv[])
 

Variables

static char * option [N_OPTIONS]
 
char * usage
 
long(* filter_funct [])(STAGE *stage, long filterNum, int64_t npoints)
 
long verbose
 

Detailed Description

General series/cascade time-domain filtering utility.

Usage

sddsdigfilter [<inputfile>] [<outputfile>] [-pipe=[input][,output]] -columns=<xcol>,<ycol> [options]

Options

-proportional=<gain> -lowpass=<gain>,<cutoff_freq> -highpass=<gain>,<cutoff_freq> -digitalfilter=<sddsfile>,<a_coeff_col>,<b_coeff_col> or -digitalfilter=[A,<a0>,<a1>,...][,B,<b0>,<b1>,...] -analogfilter=<sddsfile>,<c_coeff_col>,<d_coeff_col> or -analogfilter=[C,<c0>,<c1>,...][,D,<d0>,<d1>,...] -cascade -verbose -majorOrder=row|column

Note
Digital filter form: 
*             b0 + b1.z^-1 + ... + bn.z^-n
*     H(z) = --------------------------------
*             a0 + a1.z^-1 + ... + an.z^-n
*
Analog filter form: 
*             d0 + d1.s^1 + ... + dn.s^n
*     H(s) = ------------------------------
*             c0 + c1.s^1 + ... + cn.s^n
*
Copyright
  • (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
  • (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
License
This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.
Author
J. Carwardine, C. Saunders, R. Soliday, Xuesong Jiao, H. Shang, L. Emery

Definition in file sddsdigfilter.c.

Macro Definition Documentation

◆ FALSE

#define FALSE   0

Definition at line 51 of file sddsdigfilter.c.

◆ TRUE

#define TRUE   1

Definition at line 50 of file sddsdigfilter.c.

Enumeration Type Documentation

◆ option_type

enum option_type

Definition at line 54 of file sddsdigfilter.c.

54 {
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};

Function Documentation

◆ analog_function()

long analog_function ( STAGE * stage,
long filterNum,
int64_t npoints )

Definition at line 646 of file sddsdigfilter.c.

646 {
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}

◆ digital_function()

long digital_function ( STAGE * stage,
long filterNum,
int64_t npoints )

Definition at line 607 of file sddsdigfilter.c.

607 {
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}

◆ dpow_ri()

double dpow_ri ( double * ap,
long * bp )

Definition at line 1104 of file sddsdigfilter.c.

1104 {
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}

◆ dspbfct()

double dspbfct ( long * i1,
long * i2 )

Definition at line 1083 of file sddsdigfilter.c.

1083 {
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 */

◆ dspbiln()

void dspbiln ( double * d,
double * c,
int64_t ln,
double * b,
double * a,
double * work,
long * error )

Definition at line 819 of file sddsdigfilter.c.

819 {
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 */

◆ dspfblt()

void dspfblt ( double * d,
double * c,
long ln,
long iband,
double fln,
double fhn,
double * b,
double * a,
double * work,
long * error )

Definition at line 921 of file sddsdigfilter.c.

921 {
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 */

◆ highpass_function()

long highpass_function ( STAGE * stage,
long filterNum,
int64_t npoints )

Definition at line 536 of file sddsdigfilter.c.

536 {
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}

◆ lowpass_function()

long lowpass_function ( STAGE * stage,
long filterNum,
int64_t npoints )

Definition at line 466 of file sddsdigfilter.c.

466 {
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}

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 157 of file sddsdigfilter.c.

157 {
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}
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_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

◆ processFilterParams()

long processFilterParams ( char ** args,
long items,
char * fistStr,
char * nextStr,
ANADIG * ADptr )

Definition at line 1140 of file sddsdigfilter.c.

1140 {
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}
SDDS_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:677

◆ processFilterParamsFromFile()

long processFilterParamsFromFile ( char * filename,
char * ACColumn,
char * BDColumn,
ANADIG * ADptr )

Definition at line 1193 of file sddsdigfilter.c.

1193 {
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}

◆ proportional_function()

long proportional_function ( STAGE * stage,
long filterNum,
int64_t npoints )

Definition at line 729 of file sddsdigfilter.c.

729 {
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}

◆ response()

long response ( double * input,
double * output,
int64_t nInput,
int64_t nOutput,
double * Bcoeff,
double * Acoeff,
int64_t nBcoeffs,
int64_t nAcoeffs )

Definition at line 755 of file sddsdigfilter.c.

755 {
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}

Variable Documentation

◆ filter_funct

long(* filter_funct[])(STAGE *stage, long filterNum, int64_t npoints) ( STAGE * stage,
long filterNum,
int64_t npoints )
Initial value:
=
{
lowpass_function, highpass_function, proportional_function, analog_function, digital_function}

Definition at line 150 of file sddsdigfilter.c.

151 {
152 lowpass_function, highpass_function, proportional_function, analog_function, digital_function};

◆ option

char* option[N_OPTIONS]
static
Initial value:
= {
"lowpass",
"highpass",
"proportional",
"analogfilter",
"digitalfilter",
"cascade",
"columns",
"pipe",
"verbose",
"majorOrder"
}

Definition at line 68 of file sddsdigfilter.c.

68 {
69 "lowpass",
70 "highpass",
71 "proportional",
72 "analogfilter",
73 "digitalfilter",
74 "cascade",
75 "columns",
76 "pipe",
77 "verbose",
78 "majorOrder"
79};

◆ usage

char* usage
Initial value:
=
"\nUsage: sddsdigfilter [<inputfile>] [<outputfile>] [-pipe=[input][,output]]\n"
" -columns=<xcol>,<ycol>[,...]\n"
" [options]:\n"
" -proportional=<gain>\n"
" -lowpass=<gain>,<cutoff_freq>\n"
" -highpass=<gain>,<cutoff_freq>\n"
" -digitalfilter=<sddsfile>,<a_coeff_col>,<b_coeff_col>\n"
" or -digitalfilter=[A,<a0>,<a1>,...][,B,<b0>,<b1>,...]\n"
" -analogfilter=<sddsfile>,<c_coeff_col>,<d_coeff_col>\n"
" or -analogfilter=[C,<c0>,<c1>,...][,D,<d0>,<d1>,...]\n"
" -cascade\n"
" -verbose\n"
" -majorOrder=row|column\n"
"\nDigital filter form:\n"
" b0 + b1.z^-1 + ... + bn.z^-n\n"
" --------------------------------\n"
" a0 + a1.z^-1 + ... + an.z^-n\n"
"\nAnalog filter form:\n"
" d0 + d1.s^1 + ... + dn.s^n\n"
" ------------------------------\n"
" c0 + c1.s^1 + ... + cn.s^n\n\n"

Definition at line 81 of file sddsdigfilter.c.

◆ verbose

long verbose

Definition at line 155 of file sddsdigfilter.c.

Generated by doxygen 1.12.0