SDDSlib
Loading...
Searching...
No Matches
sddsfdfilter.c File Reference

SDDS-format frequency-domain filter program. More...

#include "mdb.h"
#include "SDDS.h"
#include "scan.h"
#include "fftpackC.h"
#include "SDDSutils.h"
#include <ctype.h>

Go to the source code of this file.

Classes

struct  THRESHOLD_FILTER
 
struct  HILO_FILTER
 
struct  NHBP_FILTER
 
struct  FILE_FILTER
 
struct  CLIP_FILTER
 
struct  FILTER
 
struct  FILTER_STAGE
 

Macros

#define FILT_START_GIVEN   0x00000001U
 
#define FILT_END_GIVEN   0x00000002U
 
#define FILT_THRES_GIVEN   0x00000004U
 
#define FILT_CENTER_GIVEN   0x00000008U
 
#define FILT_FULLWIDTH_GIVEN   0x00000010U
 
#define FILT_FREQNAME_GIVEN   0x00000020U
 
#define FILT_REALNAME_GIVEN   0x00000040U
 
#define FILT_IMAGNAME_GIVEN   0x00000080U
 
#define FILT_MAGNAME_GIVEN   0x00000100U
 
#define FILT_FRACTHRES_GIVEN   0x00000200U
 
#define FILT_LEVEL_GIVEN   0x00000400U
 
#define FILT_FILENAME_GIVEN   0x00000800U
 
#define FILT_HIGH_GIVEN   0x00001000U
 
#define FILT_LOW_GIVEN   0x00002000U
 
#define FILT_FLATWIDTH_GIVEN   0x00004000U
 
#define FL_NEWCOLUMNS   0x00001UL
 
#define FL_DIFCOLUMNS   0x00002UL
 

Enumerations

enum  option_type {
  SET_PIPE , SET_CASCADE , SET_CLIPFREQ , SET_COLUMNS ,
  SET_THRESHOLD , SET_HIGHPASS , SET_LOWPASS , SET_NOTCH ,
  SET_BANDPASS , SET_FILTERFILE , SET_NEWCOLUMNS , SET_DIFFERENCECOLUMNS ,
  SET_EXCLUDE , SET_MAJOR_ORDER , N_OPTIONS
}
 

Functions

long applyFilters (double *outputData, double *inputData, double *timeData, int64_t rows, FILTER_STAGE *filterStage, long filterStages)
 
long applyFilterStage (double *outputRI, double *inputRI, int64_t frequencies, double dfrequency, FILTER_STAGE *filterStage)
 
void addClipFilterOutput (double *outputRI, double *inputRI, int64_t frequencies, double dfrequency, CLIP_FILTER *filter)
 
void addThresholdFilterOutput (double *outputRI, double *inputRI, int64_t frequencies, double dfrequency, THRESHOLD_FILTER *filter)
 
void addHighPassFilterOutput (double *outputRI, double *inputRI, int64_t frequencies, double dfrequency, HILO_FILTER *filter)
 
void addLowPassFilterOutput (double *outputRI, double *inputRI, int64_t frequencies, double dfrequency, HILO_FILTER *filter)
 
void addNotchFilterOutput (double *outputRI, double *inputRI, int64_t frequencies, double dfrequency, NHBP_FILTER *filter)
 
void addBandPassFilterOutput (double *outputRI, double *inputRI, int64_t frequencies, double dfrequency, NHBP_FILTER *filter)
 
void addFileFilterOutput (double *outputRI, double *inputRI, int64_t frequencies, double dfrequency, FILE_FILTER *filter)
 
int64_t computeIndex (double value, double dfrequency, int64_t frequencies)
 
void addFilter (FILTER_STAGE *filterStage, long optionCode, SCANNED_ARG *scanned)
 
int main (int argc, char **argv)
 

Variables

char * option [N_OPTIONS]
 
static char * USAGE
 

Detailed Description

SDDS-format frequency-domain filter program.

This program applies various frequency-domain filters to SDDS-formatted data files. It supports multiple filtering options, including high-pass, low-pass, notch, bandpass, threshold, and file-based filters. Users can specify input and output files or use pipes.

Usage

sddsfdfilter [<inputfile>] [<outputfile>]
[-pipe=[input][,output]]
[-columns=<indep-variable>[,<depen-quantity>[,...]]]
[-exclude=<depen-quantity>[,...]]
[-clipFrequencies=[high=<number>][,low=<number>]]
[-threshold=level=<value>[,fractional][,start=<freq>][,end=<freq>]]
[-highpass=start=<freq>,end=<freq>]
[-lowpass=start=<freq>,end=<freq>]
[-notch=center=<center>,flatWidth=<value>,fullWidth=<value>]
[-bandpass=center=<center>,flatWidth=<value>,fullWidth=<value>]
[-filterFile=filename=<filename>,frequency=<columnName>{,real=<cName>,imaginary=<cName>|magnitude=<cName>}]
[-cascade]
[-newColumns]
[-differenceColumns]
[-majorOrder=row|column]
License
This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.
Author
M. Borland, R. Soliday, Xuesong Jiao, L. Emery, H. Shang

Definition in file sddsfdfilter.c.

Macro Definition Documentation

◆ FILT_CENTER_GIVEN

#define FILT_CENTER_GIVEN   0x00000008U

Definition at line 104 of file sddsfdfilter.c.

◆ FILT_END_GIVEN

#define FILT_END_GIVEN   0x00000002U

Definition at line 102 of file sddsfdfilter.c.

◆ FILT_FILENAME_GIVEN

#define FILT_FILENAME_GIVEN   0x00000800U

Definition at line 112 of file sddsfdfilter.c.

◆ FILT_FLATWIDTH_GIVEN

#define FILT_FLATWIDTH_GIVEN   0x00004000U

Definition at line 115 of file sddsfdfilter.c.

◆ FILT_FRACTHRES_GIVEN

#define FILT_FRACTHRES_GIVEN   0x00000200U

Definition at line 110 of file sddsfdfilter.c.

◆ FILT_FREQNAME_GIVEN

#define FILT_FREQNAME_GIVEN   0x00000020U

Definition at line 106 of file sddsfdfilter.c.

◆ FILT_FULLWIDTH_GIVEN

#define FILT_FULLWIDTH_GIVEN   0x00000010U

Definition at line 105 of file sddsfdfilter.c.

◆ FILT_HIGH_GIVEN

#define FILT_HIGH_GIVEN   0x00001000U

Definition at line 113 of file sddsfdfilter.c.

◆ FILT_IMAGNAME_GIVEN

#define FILT_IMAGNAME_GIVEN   0x00000080U

Definition at line 108 of file sddsfdfilter.c.

◆ FILT_LEVEL_GIVEN

#define FILT_LEVEL_GIVEN   0x00000400U

Definition at line 111 of file sddsfdfilter.c.

◆ FILT_LOW_GIVEN

#define FILT_LOW_GIVEN   0x00002000U

Definition at line 114 of file sddsfdfilter.c.

◆ FILT_MAGNAME_GIVEN

#define FILT_MAGNAME_GIVEN   0x00000100U

Definition at line 109 of file sddsfdfilter.c.

◆ FILT_REALNAME_GIVEN

#define FILT_REALNAME_GIVEN   0x00000040U

Definition at line 107 of file sddsfdfilter.c.

◆ FILT_START_GIVEN

#define FILT_START_GIVEN   0x00000001U

Definition at line 101 of file sddsfdfilter.c.

◆ FILT_THRES_GIVEN

#define FILT_THRES_GIVEN   0x00000004U

Definition at line 103 of file sddsfdfilter.c.

◆ FL_DIFCOLUMNS

#define FL_DIFCOLUMNS   0x00002UL

Definition at line 155 of file sddsfdfilter.c.

◆ FL_NEWCOLUMNS

#define FL_NEWCOLUMNS   0x00001UL

Definition at line 154 of file sddsfdfilter.c.

Enumeration Type Documentation

◆ option_type

enum option_type

Definition at line 47 of file sddsfdfilter.c.

47 {
48 SET_PIPE,
49 SET_CASCADE,
50 SET_CLIPFREQ,
51 SET_COLUMNS,
52 SET_THRESHOLD,
53 SET_HIGHPASS,
54 SET_LOWPASS,
55 SET_NOTCH,
56 SET_BANDPASS,
57 SET_FILTERFILE,
58 SET_NEWCOLUMNS,
59 SET_DIFFERENCECOLUMNS,
60 SET_EXCLUDE,
61 SET_MAJOR_ORDER,
62 N_OPTIONS
63};

Function Documentation

◆ addBandPassFilterOutput()

void addBandPassFilterOutput ( double * outputRI,
double * inputRI,
int64_t frequencies,
double dfrequency,
NHBP_FILTER * filter )

Definition at line 619 of file sddsfdfilter.c.

619 {
620 int64_t i1, i2, i;
621 double fraction, dfraction;
622
623 i1 = computeIndex(filter->center - filter->fullwidth / 2, dfrequency, frequencies);
624 i2 = computeIndex(filter->center - filter->flatwidth / 2, dfrequency, frequencies);
625 dfraction = i1 == i2 ? 0 : 1.0 / (i2 - i1);
626 fraction = 0;
627 for (i = i1; i <= i2; i++) {
628 outputRI[2 * i] += inputRI[2 * i] * fraction;
629 outputRI[2 * i + 1] += inputRI[2 * i + 1] * fraction;
630 if ((fraction += dfraction) > 1)
631 fraction = 1;
632 }
633
634 i1 = computeIndex(filter->center + filter->flatwidth / 2, dfrequency, frequencies);
635 i2 = computeIndex(filter->center + filter->fullwidth / 2, dfrequency, frequencies);
636 dfraction = i1 == i2 ? 0 : 1.0 / (i2 - i1);
637 fraction = 1;
638 for (i = i1; i <= i2; i++) {
639 outputRI[2 * i] += inputRI[2 * i] * fraction;
640 outputRI[2 * i + 1] += inputRI[2 * i + 1] * fraction;
641 if ((fraction -= dfraction) < 0)
642 fraction = 0;
643 }
644
645 i1 = computeIndex(filter->center - filter->flatwidth / 2, dfrequency, frequencies);
646 i2 = computeIndex(filter->center + filter->flatwidth / 2, dfrequency, frequencies);
647 for (i = i1; i <= i2; i++) {
648 outputRI[2 * i] += inputRI[2 * i];
649 outputRI[2 * i + 1] += inputRI[2 * i + 1];
650 }
651}

◆ addClipFilterOutput()

void addClipFilterOutput ( double * outputRI,
double * inputRI,
int64_t frequencies,
double dfrequency,
CLIP_FILTER * filter )

Definition at line 492 of file sddsfdfilter.c.

492 {
493 int64_t i1, i2;
494
495 i1 = 0;
496 i2 = frequencies - 1;
497 if (filter->flags & FILT_HIGH_GIVEN)
498 if ((i2 = frequencies - filter->high) < 0)
499 i2 = 0;
500 if (filter->flags & FILT_LOW_GIVEN)
501 if ((i1 = filter->low) >= frequencies)
502 i1 = frequencies - 1;
503 for (; i1 <= i2; i1++) {
504 outputRI[2 * i1] += inputRI[2 * i1];
505 outputRI[2 * i1 + 1] += inputRI[2 * i1 + 1];
506 }
507}

◆ addFileFilterOutput()

void addFileFilterOutput ( double * outputRI,
double * inputRI,
int64_t frequencies,
double dfrequency,
FILE_FILTER * filter )

Definition at line 653 of file sddsfdfilter.c.

653 {
654 long code;
655 int64_t i;
656 double factor, ifactor, rfactor, rdata, idata, f;
657
658 if (!filter->freqData) {
659 SDDS_DATASET SDDSin;
660 long readCode = 0;
661
662 if (!SDDS_InitializeInput(&SDDSin, filter->file) || (readCode = SDDS_ReadPage(&SDDSin)) == 0)
663 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
664 if (readCode < 0) {
665 fprintf(stderr, "error: unable to read filter file %s (sddsfdfilter)\n", filter->file);
666 exit(EXIT_FAILURE);
667 }
668 if ((filter->points = SDDS_CountRowsOfInterest(&SDDSin)) <= 0) {
669 fprintf(stderr, "error: file %s has no data on first page (sddsfdfilter)\n", filter->file);
670 exit(EXIT_FAILURE);
671 }
672 if (!(filter->freqData = SDDS_GetColumnInDoubles(&SDDSin, filter->freqName)) ||
673 (filter->magName && !(filter->magData = SDDS_GetColumnInDoubles(&SDDSin, filter->magName))) ||
674 (filter->imagName && !(filter->imagData = SDDS_GetColumnInDoubles(&SDDSin, filter->imagName))) ||
675 (filter->realName && !(filter->realData = SDDS_GetColumnInDoubles(&SDDSin, filter->realName))) ||
676 !SDDS_Terminate(&SDDSin))
677 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
678 for (i = 1; i < filter->points; i++)
679 if (filter->freqData[i - 1] >= filter->freqData[i]) {
680 fprintf(stderr, "error: frequency data not monotonically increasing for %s (sddsfdfilter)\n", filter->file);
681 exit(EXIT_FAILURE);
682 }
683 }
684 for (i = 0; i < frequencies; i++) {
685 f = i * dfrequency;
686 if (filter->freqData[0] > f || filter->freqData[filter->points - 1] < f)
687 continue;
688 if (filter->magName) {
689 factor = interp(filter->magData, filter->freqData, filter->points, f, 0, 1, &code);
690 outputRI[2 * i] += factor * inputRI[2 * i];
691 outputRI[2 * i + 1] += factor * inputRI[2 * i + 1];
692 } else {
693 rfactor = interp(filter->realData, filter->freqData, filter->points, f, 0, 1, &code);
694 ifactor = interp(filter->imagData, filter->freqData, filter->points, f, 0, 1, &code);
695 rdata = inputRI[2 * i];
696 idata = inputRI[2 * i + 1];
697 outputRI[2 * i] += rdata * rfactor - idata * ifactor;
698 outputRI[2 * i + 1] += rdata * ifactor + idata * rfactor;
699 }
700 }
701}
int64_t SDDS_CountRowsOfInterest(SDDS_DATASET *SDDS_dataset)
Counts the number of rows marked as "of interest" in the current data table.
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...
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:49
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:432
double interp(double *f, double *x, long n, double xo, long warnings, long order, long *returnCode)
Performs simple linear interpolation of data.
Definition interp.c:34

◆ addFilter()

void addFilter ( FILTER_STAGE * filterStage,
long optionCode,
SCANNED_ARG * scanned )

Definition at line 714 of file sddsfdfilter.c.

714 {
715 THRESHOLD_FILTER *thresholdFilter;
716 HILO_FILTER *hiloFilter;
717 NHBP_FILTER *nhbpFilter;
718 FILE_FILTER *fileFilter;
719 CLIP_FILTER *clipFilter;
720 FILTER *filter;
721 long items, filters;
722 char **item;
723
724 if (!(filterStage->filter = SDDS_Realloc(filterStage->filter, sizeof(*filterStage->filter) * (filterStage->filters + 1))))
725 SDDS_Bomb("allocation failure adding filter slot");
726 filter = filterStage->filter;
727 filters = filterStage->filters;
728 filterStage->filters++;
729
730 items = scanned->n_items - 1;
731 item = scanned->list + 1;
732 switch (filter[filters].filterType = optionCode) {
733 case SET_THRESHOLD:
734 if (!(thresholdFilter = filter[filters].filter = calloc(1, sizeof(THRESHOLD_FILTER))))
735 SDDS_Bomb("allocation failure allocating threshold filter");
736 if (items < 1)
737 SDDS_Bomb("invalid -threshold syntax");
738 if (!scanItemList(&thresholdFilter->flags, item, &items, 0,
739 "level", SDDS_DOUBLE, &thresholdFilter->level, 1, FILT_LEVEL_GIVEN,
740 "fractional", -1, NULL, 0, FILT_FRACTHRES_GIVEN,
741 "start", SDDS_DOUBLE, &thresholdFilter->start, 1, FILT_START_GIVEN,
742 "end", SDDS_DOUBLE, &thresholdFilter->end, 1, FILT_END_GIVEN, NULL))
743 SDDS_Bomb("invalid -threshold syntax/values");
744 if (!(thresholdFilter->flags & FILT_LEVEL_GIVEN))
745 SDDS_Bomb("supply level=<value> qualifier for -threshold");
746 if ((thresholdFilter->flags & FILT_START_GIVEN && thresholdFilter->flags & FILT_END_GIVEN) && thresholdFilter->start > thresholdFilter->end)
747 SDDS_Bomb("start > end for -threshold filter");
748 break;
749 case SET_HIGHPASS:
750 case SET_LOWPASS:
751 if (!(hiloFilter = filter[filters].filter = calloc(1, sizeof(HILO_FILTER))))
752 SDDS_Bomb("allocation failure allocating low-pass or high-pass filter");
753 if (!scanItemList(&hiloFilter->flags, item, &items, 0,
754 "start", SDDS_DOUBLE, &hiloFilter->start, 1, FILT_START_GIVEN,
755 "end", SDDS_DOUBLE, &hiloFilter->end, 1, FILT_END_GIVEN, NULL))
756 SDDS_Bomb("invalid -highpass or -lowpass syntax");
757 if (!(hiloFilter->flags & FILT_START_GIVEN) || !(hiloFilter->flags & FILT_END_GIVEN))
758 SDDS_Bomb("supply start=<value> and end=<value> qualifiers with -highpass and -lowpass");
759 break;
760 case SET_NOTCH:
761 case SET_BANDPASS:
762 if (!(nhbpFilter = filter[filters].filter = calloc(1, sizeof(NHBP_FILTER))))
763 SDDS_Bomb("allocation failure allocating notch or bandpass filter");
764 if (!scanItemList(&nhbpFilter->flags, item, &items, 0,
765 "center", SDDS_DOUBLE, &nhbpFilter->center, 1, FILT_CENTER_GIVEN,
766 "fullwidth", SDDS_DOUBLE, &nhbpFilter->fullwidth, 1, FILT_FULLWIDTH_GIVEN,
767 "flatwidth", SDDS_DOUBLE, &nhbpFilter->flatwidth, 1, FILT_FLATWIDTH_GIVEN, NULL))
768 SDDS_Bomb("invalid -notch or -bandpass syntax");
769 if (!(nhbpFilter->flags & FILT_CENTER_GIVEN) || !(nhbpFilter->flags & FILT_FLATWIDTH_GIVEN))
770 SDDS_Bomb("supply center=<value> and flatWidth=<value> qualifiers with -notch and -bandpass");
771 if (!(nhbpFilter->flags & FILT_FULLWIDTH_GIVEN))
772 nhbpFilter->fullwidth = nhbpFilter->flatwidth;
773 if (nhbpFilter->fullwidth < nhbpFilter->flatwidth)
774 SDDS_Bomb("full width may not be less than flat width for notch/bandpass filter");
775 break;
776 case SET_FILTERFILE:
777 if (!(fileFilter = filter[filters].filter = calloc(1, sizeof(FILE_FILTER))))
778 SDDS_Bomb("allocation failure allocating file filter");
779 if (!scanItemList(&fileFilter->flags, item, &items, 0,
780 "filename", SDDS_STRING, &fileFilter->file, 1, FILT_FILENAME_GIVEN,
781 "frequency", SDDS_STRING, &fileFilter->freqName, 1, FILT_FREQNAME_GIVEN,
782 "real", SDDS_STRING, &fileFilter->realName, 1, FILT_REALNAME_GIVEN,
783 "imaginary", SDDS_STRING, &fileFilter->imagName, 1, FILT_IMAGNAME_GIVEN,
784 "magnitude", SDDS_STRING, &fileFilter->magName, 1, FILT_MAGNAME_GIVEN, NULL))
785 SDDS_Bomb("invalid -filterFile syntax");
786 if (!(fileFilter->flags & FILT_FILENAME_GIVEN))
787 SDDS_Bomb("supply filename=<string> with -filterFile");
788 if (!(fileFilter->flags & FILT_FREQNAME_GIVEN))
789 SDDS_Bomb("supply frequency=<columnName> with -filterFile");
790 if (!(fileFilter->flags & FILT_MAGNAME_GIVEN) && !(fileFilter->flags & FILT_REALNAME_GIVEN && fileFilter->flags & FILT_IMAGNAME_GIVEN))
791 SDDS_Bomb("supply either magnitude=<columnName>, or real=<columnName> and imag=<columnName>, with -filterFile");
792 if (fileFilter->flags & FILT_MAGNAME_GIVEN && (fileFilter->flags & FILT_REALNAME_GIVEN || fileFilter->flags & FILT_IMAGNAME_GIVEN))
793 SDDS_Bomb("magnitude=<columnName> is incompatible with real=<columnName> and imag=<columnName> for -filterFile");
794 fileFilter->freqData = NULL; /* Indicates no data present */
795 break;
796 case SET_CLIPFREQ:
797 if (!(clipFilter = filter[filters].filter = calloc(1, sizeof(CLIP_FILTER))))
798 SDDS_Bomb("allocation failure allocating clip filter");
799 if (!scanItemList(&clipFilter->flags, item, &items, 0,
800 "high", SDDS_LONG, &clipFilter->high, 1, FILT_HIGH_GIVEN,
801 "low", SDDS_LONG, &clipFilter->low, 1, FILT_LOW_GIVEN, NULL))
802 SDDS_Bomb("invalid -clipFrequencies syntax");
803 if (!(clipFilter->flags & FILT_HIGH_GIVEN) && !(clipFilter->flags & FILT_LOW_GIVEN))
804 SDDS_Bomb("supply at least one of high=<number> or low=<number> with -clipFrequencies");
805 break;
806 default:
807 SDDS_Bomb("invalid filter code---this shouldn't happen");
808 break;
809 }
810}
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:342
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:677
#define SDDS_STRING
Identifier for the string data type.
Definition SDDStypes.h:85
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
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.

◆ addHighPassFilterOutput()

void addHighPassFilterOutput ( double * outputRI,
double * inputRI,
int64_t frequencies,
double dfrequency,
HILO_FILTER * filter )

Definition at line 542 of file sddsfdfilter.c.

542 {
543 int64_t i, i1, i2;
544 double fraction, dfraction;
545
546 i1 = computeIndex(filter->start, dfrequency, frequencies);
547 i2 = computeIndex(filter->end, dfrequency, frequencies);
548
549 dfraction = i1 == i2 ? 0 : 1.0 / (i2 - i1);
550 fraction = 0;
551 for (i = i1; i <= i2; i++) {
552 outputRI[2 * i] += inputRI[2 * i] * fraction;
553 outputRI[2 * i + 1] += inputRI[2 * i + 1] * fraction;
554 if ((fraction += dfraction) > 1)
555 fraction = 1;
556 }
557 for (; i < frequencies; i++) {
558 outputRI[2 * i] += inputRI[2 * i];
559 outputRI[2 * i + 1] += inputRI[2 * i + 1];
560 }
561}

◆ addLowPassFilterOutput()

void addLowPassFilterOutput ( double * outputRI,
double * inputRI,
int64_t frequencies,
double dfrequency,
HILO_FILTER * filter )

Definition at line 563 of file sddsfdfilter.c.

563 {
564 int64_t i, i1, i2;
565 double fraction, dfraction;
566
567 i1 = computeIndex(filter->start, dfrequency, frequencies);
568 i2 = computeIndex(filter->end, dfrequency, frequencies);
569
570 dfraction = i1 == i2 ? 0 : 1.0 / (i2 - i1);
571 fraction = 1;
572 for (i = i1; i <= i2; i++) {
573 outputRI[2 * i] += inputRI[2 * i] * fraction;
574 outputRI[2 * i + 1] += inputRI[2 * i + 1] * fraction;
575 if ((fraction -= dfraction) < 0)
576 fraction = 0;
577 }
578 for (i = 0; i < i1; i++) {
579 outputRI[2 * i] += inputRI[2 * i];
580 outputRI[2 * i + 1] += inputRI[2 * i + 1];
581 }
582}

◆ addNotchFilterOutput()

void addNotchFilterOutput ( double * outputRI,
double * inputRI,
int64_t frequencies,
double dfrequency,
NHBP_FILTER * filter )

Definition at line 584 of file sddsfdfilter.c.

584 {
585 int64_t i1, i2, i;
586 double fraction, dfraction;
587
588 i1 = computeIndex(filter->center - filter->fullwidth / 2, dfrequency, frequencies);
589 i2 = computeIndex(filter->center - filter->flatwidth / 2, dfrequency, frequencies);
590 for (i = 0; i < i1; i++) {
591 outputRI[2 * i] += inputRI[2 * i];
592 outputRI[2 * i + 1] += inputRI[2 * i + 1];
593 }
594 dfraction = i1 == i2 ? 0 : 1.0 / (i2 - i1);
595 fraction = 1;
596 for (i = i1; i <= i2; i++) {
597 outputRI[2 * i] += inputRI[2 * i] * fraction;
598 outputRI[2 * i + 1] += inputRI[2 * i + 1] * fraction;
599 if ((fraction -= dfraction) < 0)
600 fraction = 0;
601 }
602
603 i1 = computeIndex(filter->center + filter->flatwidth / 2, dfrequency, frequencies);
604 i2 = computeIndex(filter->center + filter->fullwidth / 2, dfrequency, frequencies);
605 for (i = i2; i < frequencies; i++) {
606 outputRI[2 * i] += inputRI[2 * i];
607 outputRI[2 * i + 1] += inputRI[2 * i + 1];
608 }
609 dfraction = i1 == i2 ? 0 : 1.0 / (i2 - i1);
610 fraction = 0;
611 for (i = i1; i <= i2; i++) {
612 outputRI[2 * i] += inputRI[2 * i] * fraction;
613 outputRI[2 * i + 1] += inputRI[2 * i + 1] * fraction;
614 if ((fraction += dfraction) > 1)
615 fraction = 1;
616 }
617}

◆ addThresholdFilterOutput()

void addThresholdFilterOutput ( double * outputRI,
double * inputRI,
int64_t frequencies,
double dfrequency,
THRESHOLD_FILTER * filter )

Definition at line 509 of file sddsfdfilter.c.

509 {
510 int64_t i, i1, i2;
511 double level2, level2q, level2o;
512 double max2, mag2;
513
514 i1 = 0;
515 if (filter->flags & FILT_START_GIVEN)
516 i1 = computeIndex(filter->start, dfrequency, frequencies);
517 i2 = frequencies - 1;
518 if (filter->flags & FILT_END_GIVEN)
519 i2 = computeIndex(filter->end, dfrequency, frequencies);
520
521 if (filter->flags & FILT_FRACTHRES_GIVEN) {
522 max2 = -DBL_MAX;
523 for (i = i1; i <= i2; i++)
524 if ((mag2 = inputRI[2 * i] * inputRI[2 * i] + inputRI[2 * i + 1] * inputRI[2 * i + 1]) > max2)
525 max2 = mag2;
526 level2o = max2 * sqr(filter->level);
527 } else
528 level2o = sqr(filter->level);
529 level2q = level2o / 4;
530
531 for (i = i1; i <= i2; i++) {
532 level2 = level2q; /* only half the amplitude is in positive frequencies */
533 if (i == 0 || i == frequencies - 1)
534 level2 = level2o;
535 if ((inputRI[2 * i] * inputRI[2 * i] + inputRI[2 * i + 1] * inputRI[2 * i + 1]) >= level2) {
536 outputRI[2 * i] += inputRI[2 * i];
537 outputRI[2 * i + 1] += inputRI[2 * i + 1];
538 }
539 }
540}

◆ applyFilters()

long applyFilters ( double * outputData,
double * inputData,
double * timeData,
int64_t rows,
FILTER_STAGE * filterStage,
long filterStages )

Definition at line 410 of file sddsfdfilter.c.

410 {
411 long stage;
412 int64_t frequencies, row;
413 double *realimagInput, *realimagOutput;
414 double length, dfrequency, factor;
415 realimagOutput = NULL;
416 if (!(realimagInput = (double *)malloc(sizeof(*realimagInput) * (rows + 2))) ||
417 !(realimagOutput = (double *)malloc(sizeof(*realimagOutput) * (rows + 2))))
418 SDDS_Bomb("allocation failure");
419
420 /* Input data is real and has "rows" rows */
421 /* Result is interleaved real and imaginary */
422 realFFT2(realimagInput, inputData, rows, 0);
423 frequencies = rows / 2 + 1;
424
425 /* Length is the assumed length of the periodic signal,
426 which is one interval longer than the actual data entered */
427 length = ((double)rows) * (timeData[rows - 1] - timeData[0]) / ((double)rows - 1.0);
428 dfrequency = factor = 1.0 / length;
429
430 for (stage = 0; stage < filterStages; stage++) {
431 if (!applyFilterStage(realimagOutput, realimagInput, frequencies, dfrequency, filterStage + stage))
432 return 0;
433 SWAP_PTR(realimagOutput, realimagInput);
434 }
435
436 /* Input is still interleaved real and imaginary. The flag
437 INVERSE_FFT ensures that the array is interpreted correctly. */
438 /* Output is interleaved real and imaginary */
439#if DEBUG
440 for (row = 0; row < rows; row++) {
441 fprintf(stdout, "Real: %f, Imag: %f\n", realimagInput[2 * row], realimagInput[2 * row + 1]);
442 }
443#endif
444 realFFT2(realimagOutput, realimagInput, rows, INVERSE_FFT);
445
446 for (row = 0; row < rows; row++)
447 outputData[row] = realimagOutput[row];
448
449 free(realimagOutput);
450 free(realimagInput);
451
452 return 1;
453}

◆ applyFilterStage()

long applyFilterStage ( double * outputRI,
double * inputRI,
int64_t frequencies,
double dfrequency,
FILTER_STAGE * filterStage )

Definition at line 455 of file sddsfdfilter.c.

455 {
456 int64_t i, ifilter;
457
458 for (i = 0; i < 2 * frequencies; i++)
459 outputRI[i] = 0;
460
461 for (ifilter = 0; ifilter < filterStage->filters; ifilter++) {
462 switch (filterStage->filter[ifilter].filterType) {
463 case SET_CLIPFREQ:
464 addClipFilterOutput(outputRI, inputRI, frequencies, dfrequency, filterStage->filter[ifilter].filter);
465 break;
466 case SET_THRESHOLD:
467 addThresholdFilterOutput(outputRI, inputRI, frequencies, dfrequency, filterStage->filter[ifilter].filter);
468 break;
469 case SET_HIGHPASS:
470 addHighPassFilterOutput(outputRI, inputRI, frequencies, dfrequency, filterStage->filter[ifilter].filter);
471 break;
472 case SET_LOWPASS:
473 addLowPassFilterOutput(outputRI, inputRI, frequencies, dfrequency, filterStage->filter[ifilter].filter);
474 break;
475 case SET_NOTCH:
476 addNotchFilterOutput(outputRI, inputRI, frequencies, dfrequency, filterStage->filter[ifilter].filter);
477 break;
478 case SET_BANDPASS:
479 addBandPassFilterOutput(outputRI, inputRI, frequencies, dfrequency, filterStage->filter[ifilter].filter);
480 break;
481 case SET_FILTERFILE:
482 addFileFilterOutput(outputRI, inputRI, frequencies, dfrequency, filterStage->filter[ifilter].filter);
483 break;
484 default:
485 SDDS_Bomb("unknown filter type in applyFilterStage--this shouldn't happen");
486 break;
487 }
488 }
489 return 1;
490}

◆ computeIndex()

int64_t computeIndex ( double value,
double dfrequency,
int64_t frequencies )

Definition at line 703 of file sddsfdfilter.c.

703 {
704 int64_t i1;
705
706 i1 = value / dfrequency + 0.5;
707 if (i1 < 0)
708 i1 = 0;
709 if (i1 > frequencies - 1)
710 i1 = frequencies - 1;
711 return i1;
712}

◆ main()

int main ( int argc,
char ** argv )

Definition at line 169 of file sddsfdfilter.c.

169 {
170 int iArg;
171 char **outputColumn, **difColumn;
172 char *indepColumn, **depenColumn, **exclude;
173 long depenColumns, excludes;
174 char *input, *output;
175 long i, readCode, optionCode;
176 int64_t rows;
177 unsigned long flags, pipeFlags, majorOrderFlag;
178 SCANNED_ARG *scanned;
179 SDDS_DATASET SDDSin, SDDSout;
180 double *timeData, *inputData, *outputData;
181 FILTER_STAGE *filterStage;
182 long filterStages, totalFilters;
183 short columnMajorOrder = -1;
184
186 argc = scanargs(&scanned, argc, argv);
187 if (argc < 3 || argc > (3 + N_OPTIONS))
188 bomb(NULL, USAGE);
189
190 output = input = NULL;
191 flags = pipeFlags = 0;
192 indepColumn = NULL;
193 depenColumn = exclude = NULL;
194 depenColumns = excludes = 0;
195 if (!(filterStage = (FILTER_STAGE *)calloc(1, sizeof(*filterStage))))
196 SDDS_Bomb("allocation failure");
197
198 filterStage->filter = NULL;
199 filterStage->filters = 0;
200 filterStages = 1;
201 totalFilters = 0;
202
203 for (iArg = 1; iArg < argc; iArg++) {
204 if (scanned[iArg].arg_type == OPTION) {
205 /* Process options here */
206 switch (optionCode = match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
207 case SET_MAJOR_ORDER:
208 majorOrderFlag = 0;
209 scanned[iArg].n_items--;
210 if (scanned[iArg].n_items > 0 && (!scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
211 "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
212 "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
213 SDDS_Bomb("invalid -majorOrder syntax/values");
214 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
215 columnMajorOrder = 1;
216 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
217 columnMajorOrder = 0;
218 break;
219 case SET_PIPE:
220 if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
221 SDDS_Bomb("invalid -pipe syntax");
222 break;
223 case SET_COLUMNS:
224 if (indepColumn)
225 SDDS_Bomb("only one -columns option may be given");
226 if (scanned[iArg].n_items < 2)
227 SDDS_Bomb("invalid -columns syntax");
228 indepColumn = scanned[iArg].list[1];
229 if (scanned[iArg].n_items >= 2) {
230 depenColumn = tmalloc(sizeof(*depenColumn) * (depenColumns = scanned[iArg].n_items - 2));
231 for (i = 0; i < depenColumns; i++)
232 depenColumn[i] = scanned[iArg].list[i + 2];
233 }
234 break;
235 case SET_THRESHOLD:
236 case SET_HIGHPASS:
237 case SET_LOWPASS:
238 case SET_NOTCH:
239 case SET_BANDPASS:
240 case SET_FILTERFILE:
241 case SET_CLIPFREQ:
242 addFilter(filterStage + filterStages - 1, optionCode, scanned + iArg);
243 totalFilters++;
244 break;
245 case SET_CASCADE:
246 if (filterStages == 0)
247 SDDS_Bomb("-cascade option precedes all filter definitions");
248 if (!(filterStage = SDDS_Realloc(filterStage, (filterStages + 1) * sizeof(*filterStage))))
249 SDDS_Bomb("allocation failure");
250 filterStage[filterStages].filter = NULL;
251 filterStage[filterStages].filters = 0;
252 filterStages++;
253 break;
254 case SET_NEWCOLUMNS:
255 flags |= FL_NEWCOLUMNS;
256 break;
257 case SET_DIFFERENCECOLUMNS:
258 flags |= FL_DIFCOLUMNS;
259 break;
260 case SET_EXCLUDE:
261 if (scanned[iArg].n_items < 2)
262 SDDS_Bomb("invalid -exclude syntax");
263 moveToStringArray(&exclude, &excludes, scanned[iArg].list + 1, scanned[iArg].n_items - 1);
264 break;
265 default:
266 fprintf(stderr, "error: unknown/ambiguous option: %s (%s)\n", scanned[iArg].list[0], argv[0]);
267 exit(EXIT_FAILURE);
268 break;
269 }
270 } else {
271 if (!input)
272 input = scanned[iArg].list[0];
273 else if (!output)
274 output = scanned[iArg].list[0];
275 else
276 SDDS_Bomb("too many filenames seen");
277 }
278 }
279
280 processFilenames("sddsfdfilter", &input, &output, pipeFlags, 0, NULL);
281
282 if (!totalFilters)
283 fputs("warning: no filters specified (sddsfdfilter)\n", stderr);
284
285 if (!indepColumn)
286 SDDS_Bomb("supply the independent column name with the -columns option");
287
288 if (!SDDS_InitializeInput(&SDDSin, input))
289 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
290
291 if (SDDS_CheckColumn(&SDDSin, indepColumn, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OKAY)
292 exit(EXIT_FAILURE);
293
294 excludes = appendToStringArray(&exclude, excludes, indepColumn);
295 if (!depenColumns)
296 depenColumns = appendToStringArray(&depenColumn, depenColumns, "*");
297
298 if ((depenColumns = expandColumnPairNames(&SDDSin, &depenColumn, NULL, depenColumns, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
299 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
300 SDDS_Bomb("No quantities selected to filter");
301 }
302
303 if (!SDDS_InitializeCopy(&SDDSout, &SDDSin, output, "w"))
304 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
305
306 if (columnMajorOrder != -1)
307 SDDSout.layout.data_mode.column_major = columnMajorOrder;
308 else
309 SDDSout.layout.data_mode.column_major = SDDSin.layout.data_mode.column_major;
310
311 if (flags & FL_NEWCOLUMNS) {
312 outputColumn = tmalloc(sizeof(*outputColumn) * depenColumns);
313 for (i = 0; i < depenColumns; i++) {
314 outputColumn[i] = tmalloc(sizeof(**outputColumn) * (strlen(depenColumn[i]) + 1 + strlen("Filtered")));
315 sprintf(outputColumn[i], "%sFiltered", depenColumn[i]);
316 if (!SDDS_TransferColumnDefinition(&SDDSout, &SDDSin, depenColumn[i], outputColumn[i]))
317 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
318 }
319 } else
320 outputColumn = depenColumn;
321
322 difColumn = NULL;
323 if (flags & FL_DIFCOLUMNS) {
324 difColumn = tmalloc(sizeof(*difColumn) * depenColumns);
325 for (i = 0; i < depenColumns; i++) {
326 difColumn[i] = tmalloc(sizeof(**difColumn) * (strlen(depenColumn[i]) + 1 + strlen("Difference")));
327 sprintf(difColumn[i], "%sDifference", depenColumn[i]);
328 if (!SDDS_TransferColumnDefinition(&SDDSout, &SDDSin, depenColumn[i], difColumn[i]))
329 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
330 }
331 }
332
333 if (!SDDS_WriteLayout(&SDDSout))
334 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
335
336 outputData = NULL;
337 while ((readCode = SDDS_ReadPage(&SDDSin)) > 0) {
338 if (!SDDS_CopyPage(&SDDSout, &SDDSin))
339 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
340 if ((rows = SDDS_CountRowsOfInterest(&SDDSin)) < 0)
341 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
342
343 if (rows) {
344 if (!(timeData = SDDS_GetColumnInDoubles(&SDDSin, indepColumn)))
345 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
346 if (!(outputData = SDDS_Realloc(outputData, sizeof(*outputData) * rows)))
347 SDDS_Bomb("allocation failure");
348 for (i = 0; i < depenColumns; i++) {
349 if (!(inputData = SDDS_GetColumnInDoubles(&SDDSin, depenColumn[i])))
350 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
351 if (!applyFilters(outputData, inputData, timeData, rows, filterStage, filterStages))
352 exit(EXIT_FAILURE);
353 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_BY_NAME, outputData, rows, outputColumn[i]))
354 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
355 if (flags & FL_DIFCOLUMNS) {
356 int64_t j;
357 for (j = 0; j < rows; j++)
358 outputData[j] = inputData[j] - outputData[j];
359 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_BY_NAME, outputData, rows, difColumn[i]))
360 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
361 }
362 free(inputData);
363 }
364 free(timeData);
365 }
366 if (!SDDS_WritePage(&SDDSout))
367 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
368 }
369
370 free(outputData);
371 for (i = 0; i < depenColumns; i++) {
372 free(depenColumn[i]);
373 if (flags & FL_NEWCOLUMNS)
374 free(outputColumn[i]);
375 if (flags & FL_DIFCOLUMNS)
376 free(difColumn[i]);
377 }
378 for (i = 0; i < excludes; i++)
379 free(exclude[i]);
380
381 free(indepColumn);
382 if (flags & FL_NEWCOLUMNS)
383 free(outputColumn);
384 free(depenColumn);
385 if (flags & FL_DIFCOLUMNS)
386 free(difColumn);
387 free(exclude);
388 for (i = 0; i < filterStages; i++) {
389 long j;
390 for (j = 0; j < filterStage[i].filters; j++) {
391 switch (filterStage[i].filter[j].filterType) {
392 case SET_FILTERFILE:
393 free(((FILE_FILTER *)(filterStage[i].filter[j].filter))->freqData);
394 free(((FILE_FILTER *)(filterStage[i].filter[j].filter))->magData);
395 free(((FILE_FILTER *)(filterStage[i].filter[j].filter))->imagData);
396 free(((FILE_FILTER *)(filterStage[i].filter[j].filter))->realData);
397 break;
398 default:
399 break;
400 }
401 }
402 }
403
404 if (!SDDS_Terminate(&SDDSout) || !SDDS_Terminate(&SDDSin))
405 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
406
407 return EXIT_SUCCESS;
408}
int32_t SDDS_InitializeCopy(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source, char *filename, char *filemode)
Definition SDDS_copy.c:40
int32_t SDDS_CopyPage(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
Definition SDDS_copy.c:578
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.
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
int32_t SDDS_TransferColumnDefinition(SDDS_DATASET *target, SDDS_DATASET *source, char *name, char *newName)
Transfers a column definition from a source dataset to a target dataset.
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.
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:288
#define SDDS_ANY_NUMERIC_TYPE
Special identifier used by SDDS_Check*() routines to accept any numeric type.
Definition SDDStypes.h:157
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:59
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
Definition bomb.c:26
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.
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
Definition scanargs.c:36
long processPipeOption(char **item, long items, unsigned long *flags)
Definition scanargs.c:356
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:390

Variable Documentation

◆ option

char* option[N_OPTIONS]
Initial value:
= {
"pipe",
"cascade",
"clip",
"columns",
"threshold",
"highpass",
"lowpass",
"notch",
"bandpass",
"filterfile",
"newcolumns",
"differencecolumns",
"exclude",
"majorOrder",
}

Definition at line 65 of file sddsfdfilter.c.

65 {
66 "pipe",
67 "cascade",
68 "clip",
69 "columns",
70 "threshold",
71 "highpass",
72 "lowpass",
73 "notch",
74 "bandpass",
75 "filterfile",
76 "newcolumns",
77 "differencecolumns",
78 "exclude",
79 "majorOrder",
80};

◆ USAGE

char* USAGE
static
Initial value:
=
"Usage: sddsfdfilter [<inputfile>] [<outputfile>]\n"
" [-pipe=[input][,output]]\n"
" [-columns=<indep-variable>[,<depen-quantity>[,...]]]\n"
" [-exclude=<depen-quantity>[,...]]\n"
" [-clipFrequencies=[high=<number>][,low=<number>]]\n"
" [-threshold=level=<value>[,fractional][,start=<freq>][,end=<freq>]]\n"
" [-highpass=start=<freq>,end=<freq>]\n"
" [-lowpass=start=<freq>,end=<freq>]\n"
" [-notch=center=<center>,flatWidth=<value>,fullWidth=<value>]\n"
" [-bandpass=center=<center>,flatWidth=<value>,fullWidth=<value>]\n"
" [-filterFile=filename=<filename>,frequency=<columnName>{,real=<cName>,imaginary=<cName>|magnitude=<cName>}]\n"
" [-cascade]\n"
" [-newColumns]\n"
" [-differenceColumns]\n"
" [-majorOrder=row|column]\n\n"
"Program by Michael Borland. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n"

Definition at line 83 of file sddsfdfilter.c.