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

SDDS-format FFT (Fast Fourier Transform) 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.

Macros

#define FL_TRUNCATE   0x0001
 
#define FL_PADWITHZEROES   0x0002
 
#define FL_NORMALIZE   0x0004
 
#define FL_SUPPRESSAVERAGE   0x0008
 
#define FL_FULLOUTPUT   0x0010
 
#define FL_MAKEFREQDATA   0x0020
 
#define FL_PSDOUTPUT   0x0040
 
#define FL_PSDINTEGOUTPUT   0x0080
 
#define FL_PSDRINTEGOUTPUT   0x0100
 
#define FL_FULLOUTPUT_FOLDED   0x0200
 
#define FL_FULLOUTPUT_UNFOLDED   0x0400
 
#define FL_COMPLEXINPUT_FOLDED   0x0800
 
#define FL_COMPLEXINPUT_UNFOLDED   0x1000
 
#define FL_UNWRAP_PHASE   0x2000
 
#define WINDOW_HANNING   0
 
#define WINDOW_WELCH   1
 
#define WINDOW_PARZEN   2
 
#define WINDOW_HAMMING   3
 
#define WINDOW_FLATTOP   4
 
#define WINDOW_GAUSSIAN   5
 
#define WINDOW_NONE   6
 
#define N_WINDOW_TYPES   7
 

Enumerations

enum  option_type {
  SET_WINDOW , SET_NORMALIZE , SET_PADWITHZEROES , SET_TRUNCATE ,
  SET_SUPPRESSAVERAGE , SET_SAMPLEINTERVAL , SET_COLUMNS , SET_FULLOUTPUT ,
  SET_PIPE , SET_PSDOUTPUT , SET_EXCLUDE , SET_NOWARNINGS ,
  SET_COMPLEXINPUT , SET_INVERSE , SET_MAJOR_ORDER , N_OPTIONS
}
 

Functions

int64_t greatestProductOfSmallPrimes (int64_t rows)
 
long process_data (SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, double *tdata, int64_t rows, int64_t rowsToUse, char *depenQuantity, char *depenQuantity2, unsigned long flags, long windowType, int64_t sampleInterval, long correctWindowEffects, long inverse, double rintegCutOffFreq, double unwrapLimit)
 
long create_fft_frequency_column (SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *timeName, char *freqUnits, long inverse)
 
long create_fft_columns (SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *origName, char *indepName, char *freqUnits, long full_output, unsigned long psd_output, long complexInput, long inverse, long unwrap_phase)
 
long create_fft_parameters (SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, char *indepName, char *freqUnits)
 
char * makeFrequencyUnits (SDDS_DATASET *SDDSin, char *indepName)
 
long expandComplexColumnPairNames (SDDS_DATASET *SDDSin, char **name, char ***realName, char ***imagName, long names, char **excludeName, long excludeNames, long typeMode, long typeValue)
 
int main (int argc, char **argv)
 
void moveToStringArrayComplex (char ***targetReal, char ***targetImag, long *targets, char **sourceReal, char **sourceImag, long sources)
 

Variables

char * option [N_OPTIONS]
 
char * window_type [N_WINDOW_TYPES]
 
static char * USAGE1
 
static char * USAGE2
 
static long psdOffset
 
static long argOffset
 
static long realOffset
 
static long imagOffset
 
static long fftOffset = -1
 
static long psdIntOffset
 
static long psdIntPowerOffset
 
static long unwrappedArgOffset = -1
 

Detailed Description

SDDS-format FFT (Fast Fourier Transform) program.

This program performs Fast Fourier Transforms on SDDS-formatted data files. It supports various windowing functions, normalization, padding/truncation, Power Spectral Density (PSD) outputs, and more.

Usage

sddsfft [<inputfile>] [<outputfile>]
[-pipe=[input][,output]]
[-columns=<indep-variable>[,<depen-quantity>[,...]]]
[-complexInput[=unfolded|folded]]
[-exclude=<depen-quantity>[,...]]
[-window[={hanning|welch|parzen|hamming|flattop|gaussian|none}[,correct]]]
[-sampleInterval=<number>]
[-normalize]
[-fullOutput[=unfolded|folded],unwrapLimit=<value>]
[-psdOutput[=plain][,{integrated|rintegrated[=<cutoff>]}]]
[-inverse]
[-padwithzeroes[=exponent] | -truncate]
[-suppressaverage]
[-noWarnings]
[-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, C. Saunders, R. Soliday, L. Emery

Definition in file sddsfft.c.

Macro Definition Documentation

◆ FL_COMPLEXINPUT_FOLDED

#define FL_COMPLEXINPUT_FOLDED   0x0800

Definition at line 95 of file sddsfft.c.

◆ FL_COMPLEXINPUT_UNFOLDED

#define FL_COMPLEXINPUT_UNFOLDED   0x1000

Definition at line 96 of file sddsfft.c.

◆ FL_FULLOUTPUT

#define FL_FULLOUTPUT   0x0010

Definition at line 88 of file sddsfft.c.

◆ FL_FULLOUTPUT_FOLDED

#define FL_FULLOUTPUT_FOLDED   0x0200

Definition at line 93 of file sddsfft.c.

◆ FL_FULLOUTPUT_UNFOLDED

#define FL_FULLOUTPUT_UNFOLDED   0x0400

Definition at line 94 of file sddsfft.c.

◆ FL_MAKEFREQDATA

#define FL_MAKEFREQDATA   0x0020

Definition at line 89 of file sddsfft.c.

◆ FL_NORMALIZE

#define FL_NORMALIZE   0x0004

Definition at line 86 of file sddsfft.c.

◆ FL_PADWITHZEROES

#define FL_PADWITHZEROES   0x0002

Definition at line 85 of file sddsfft.c.

◆ FL_PSDINTEGOUTPUT

#define FL_PSDINTEGOUTPUT   0x0080

Definition at line 91 of file sddsfft.c.

◆ FL_PSDOUTPUT

#define FL_PSDOUTPUT   0x0040

Definition at line 90 of file sddsfft.c.

◆ FL_PSDRINTEGOUTPUT

#define FL_PSDRINTEGOUTPUT   0x0100

Definition at line 92 of file sddsfft.c.

◆ FL_SUPPRESSAVERAGE

#define FL_SUPPRESSAVERAGE   0x0008

Definition at line 87 of file sddsfft.c.

◆ FL_TRUNCATE

#define FL_TRUNCATE   0x0001

Definition at line 84 of file sddsfft.c.

◆ FL_UNWRAP_PHASE

#define FL_UNWRAP_PHASE   0x2000

Definition at line 97 of file sddsfft.c.

◆ N_WINDOW_TYPES

#define N_WINDOW_TYPES   7

Definition at line 106 of file sddsfft.c.

◆ WINDOW_FLATTOP

#define WINDOW_FLATTOP   4

Definition at line 103 of file sddsfft.c.

◆ WINDOW_GAUSSIAN

#define WINDOW_GAUSSIAN   5

Definition at line 104 of file sddsfft.c.

◆ WINDOW_HAMMING

#define WINDOW_HAMMING   3

Definition at line 102 of file sddsfft.c.

◆ WINDOW_HANNING

#define WINDOW_HANNING   0

Definition at line 99 of file sddsfft.c.

◆ WINDOW_NONE

#define WINDOW_NONE   6

Definition at line 105 of file sddsfft.c.

◆ WINDOW_PARZEN

#define WINDOW_PARZEN   2

Definition at line 101 of file sddsfft.c.

◆ WINDOW_WELCH

#define WINDOW_WELCH   1

Definition at line 100 of file sddsfft.c.

Enumeration Type Documentation

◆ option_type

enum option_type

Definition at line 47 of file sddsfft.c.

47 {
48 SET_WINDOW,
49 SET_NORMALIZE,
50 SET_PADWITHZEROES,
51 SET_TRUNCATE,
52 SET_SUPPRESSAVERAGE,
53 SET_SAMPLEINTERVAL,
54 SET_COLUMNS,
55 SET_FULLOUTPUT,
56 SET_PIPE,
57 SET_PSDOUTPUT,
58 SET_EXCLUDE,
59 SET_NOWARNINGS,
60 SET_COMPLEXINPUT,
61 SET_INVERSE,
62 SET_MAJOR_ORDER,
63 N_OPTIONS
64};

Function Documentation

◆ create_fft_columns()

long create_fft_columns ( SDDS_DATASET * SDDSout,
SDDS_DATASET * SDDSin,
char * origName,
char * indepName,
char * freqUnits,
long full_output,
unsigned long psd_output,
long complexInput,
long inverse,
long unwrap_phase )

Definition at line 907 of file sddsfft.c.

907 {
908 char s[SDDS_MAXLINE];
909 char *origUnits, *origSymbol;
910 char *description, *name, *symbol, *units;
911 long index0, index1;
912 long offset = 0;
913
914 if (complexInput)
915 offset = 4;
916 if (SDDS_GetColumnInformation(SDDSin, "units", &origUnits, SDDS_GET_BY_NAME, origName) != SDDS_STRING ||
917 SDDS_GetColumnInformation(SDDSin, "symbol", &origSymbol, SDDS_GET_BY_NAME, origName) != SDDS_STRING)
918 return 0;
919 if (!inverse)
920 sprintf(s, "FFT%s", origName + offset);
921 else {
922 if (strncmp(origName, "FFT", 3) == 0)
923 offset = 3;
924 else if (strncmp(origName, "RealFFT", 7) == 0)
925 offset = 7;
926 else
927 offset = 0;
928 sprintf(s, "%s", origName + offset);
929 }
930 SDDS_CopyString(&name, s);
931 if (!origSymbol)
932 SDDS_CopyString(&origSymbol, origName + offset);
933 sprintf(s, "FFT %s", origSymbol);
934 SDDS_CopyString(&symbol, s);
935
936 sprintf(s, "Amplitude of FFT of %s", origSymbol);
937 SDDS_CopyString(&description, s);
938
939 if (SDDS_NumberOfErrors() || (index0 = SDDS_DefineColumn(SDDSout, name, symbol, origUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
940 return 0;
941 free(name);
942 free(symbol);
943 free(description);
944
945 if (fftOffset == -1)
946 fftOffset = 0;
947
948 if (psd_output & FL_PSDOUTPUT) {
949 if (origUnits && !SDDS_StringIsBlank(origUnits)) {
950 if (freqUnits && !SDDS_StringIsBlank(freqUnits)) {
951 sprintf(s, "(%s)$a2$n/(%s)", origUnits, freqUnits);
952 } else
953 sprintf(s, "(%s)$a2$n", origUnits);
954 SDDS_CopyString(&units, s);
955 } else
956 units = NULL;
957
958 sprintf(s, "PSD%s", origName + offset);
959 SDDS_CopyString(&name, s);
960
961 if (!origSymbol)
962 SDDS_CopyString(&origSymbol, origName + offset);
963 sprintf(s, "PSD %s", origSymbol);
964 SDDS_CopyString(&symbol, s);
965
966 sprintf(s, "PSD of %s", origSymbol);
967 SDDS_CopyString(&description, s);
968
969 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, units, description, NULL, SDDS_DOUBLE, 0)) < 0)
970 return 0;
971 psdOffset = index1 - index0;
972 free(name);
973 if (units)
974 free(units);
975 free(symbol);
976 free(description);
977 }
978
979 if (psd_output & (FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT)) {
980 if (origUnits && !SDDS_StringIsBlank(origUnits)) {
981 SDDS_CopyString(&units, origUnits);
982 } else
983 units = NULL;
984
985 sprintf(s, "SqrtIntegPSD%s", origName + offset);
986 SDDS_CopyString(&name, s);
987
988 if (!origSymbol)
989 SDDS_CopyString(&origSymbol, origName + offset);
990 sprintf(s, "Sqrt Integ PSD %s", origSymbol);
991 SDDS_CopyString(&symbol, s);
992
993 sprintf(s, "Sqrt Integ PSD of %s", origSymbol);
994 SDDS_CopyString(&description, s);
995
996 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, units, description, NULL, SDDS_DOUBLE, 0)) < 0)
997 return 0;
998 psdIntOffset = index1 - index0;
999
1000 sprintf(s, "IntegPSD%s", origName + offset);
1001 SDDS_CopyString(&name, s);
1002
1003 if (!origSymbol)
1004 SDDS_CopyString(&origSymbol, origName + offset);
1005 sprintf(s, "Integ PSD %s", origSymbol);
1006 SDDS_CopyString(&symbol, s);
1007
1008 sprintf(s, "Integ PSD of %s", origSymbol);
1009 SDDS_CopyString(&description, s);
1010
1011 if (origUnits && !SDDS_StringIsBlank(origUnits)) {
1012 sprintf(s, "%sPower", origUnits);
1013 SDDS_CopyString(&units, origUnits);
1014 } else
1015 units = NULL;
1016
1017 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, units, description, NULL, SDDS_DOUBLE, 0)) < 0)
1018 return 0;
1019 psdIntPowerOffset = index1 - index0;
1020 free(name);
1021 if (units)
1022 free(units);
1023 free(symbol);
1024 free(description);
1025 }
1026
1027 if (full_output) {
1028 if (!inverse)
1029 sprintf(s, "RealFFT%s", origName + offset);
1030 else
1031 sprintf(s, "Real%s", origName + offset);
1032 SDDS_CopyString(&name, s);
1033
1034 if (!origSymbol)
1035 SDDS_CopyString(&origSymbol, origName + offset);
1036 if (!inverse)
1037 sprintf(s, "Re[FFT %s]", origSymbol);
1038 else
1039 sprintf(s, "Re[%s]", origSymbol);
1040 SDDS_CopyString(&symbol, s);
1041
1042 if (!inverse)
1043 sprintf(s, "Real part of FFT of %s", origSymbol);
1044 else
1045 sprintf(s, "Real part of %s", origSymbol);
1046 SDDS_CopyString(&description, s);
1047
1048 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, origUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
1049 return 0;
1050 realOffset = index1 - index0;
1051 free(name);
1052 free(symbol);
1053 free(description);
1054
1055 if (!inverse)
1056 sprintf(s, "ImagFFT%s", origName + offset);
1057 else
1058 sprintf(s, "Imag%s", origName + offset);
1059 SDDS_CopyString(&name, s);
1060
1061 if (!origSymbol)
1062 SDDS_CopyString(&origSymbol, origName + offset);
1063 if (!inverse)
1064 sprintf(s, "Im[FFT %s]", origSymbol);
1065 else
1066 sprintf(s, "Im[%s]", origSymbol);
1067 SDDS_CopyString(&symbol, s);
1068
1069 if (!inverse)
1070 sprintf(s, "Imaginary part of FFT of %s", origSymbol);
1071 else
1072 sprintf(s, "Imaginary part of %s", origSymbol);
1073 SDDS_CopyString(&description, s);
1074
1075 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, origUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
1076 return 0;
1077 imagOffset = index1 - index0;
1078 free(name);
1079 free(symbol);
1080 free(description);
1081
1082 if (!inverse)
1083 sprintf(s, "ArgFFT%s", origName + offset);
1084 else
1085 sprintf(s, "Arg%s", origName + offset);
1086 SDDS_CopyString(&name, s);
1087
1088 if (!origSymbol)
1089 SDDS_CopyString(&origSymbol, origName + offset);
1090 if (!inverse)
1091 sprintf(s, "Arg[FFT %s]", origSymbol);
1092 else
1093 sprintf(s, "Arg[%s]", origSymbol);
1094 SDDS_CopyString(&symbol, s);
1095
1096 if (!inverse)
1097 sprintf(s, "Phase of FFT of %s", origSymbol);
1098 else
1099 sprintf(s, "Phase of %s", origSymbol);
1100 SDDS_CopyString(&description, s);
1101
1102 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, "degrees", description, NULL, SDDS_DOUBLE, 0)) < 0)
1103 return 0;
1104 argOffset = index1 - index0;
1105 free(name);
1106 free(symbol);
1107 free(description);
1108 if (unwrap_phase) {
1109 if (!inverse)
1110 sprintf(s, "UnwrapArgFFT%s", origName + offset);
1111 else
1112 sprintf(s, "UnwrapArg%s", origName + offset);
1113 SDDS_CopyString(&name, s);
1114
1115 if (!origSymbol)
1116 SDDS_CopyString(&origSymbol, origName + offset);
1117 if (!inverse)
1118 sprintf(s, "UnwrapArg[FFT %s]", origSymbol);
1119 else
1120 sprintf(s, "UnwrapArg[%s]", origSymbol);
1121 SDDS_CopyString(&symbol, s);
1122
1123 if (!inverse)
1124 sprintf(s, "Unwrapped Phase of FFT of %s", origSymbol);
1125 else
1126 sprintf(s, "Unwrapped Phase of %s", origSymbol);
1127 SDDS_CopyString(&description, s);
1128
1129 if (SDDS_NumberOfErrors() || (index1 = SDDS_DefineColumn(SDDSout, name, symbol, "degrees", description, NULL, SDDS_DOUBLE, 0)) < 0)
1130 return 0;
1131 unwrappedArgOffset = index1 - index0;
1132 free(name);
1133 free(symbol);
1134 free(description);
1135 }
1136 }
1137
1138 return 1;
1139}
int32_t SDDS_GetColumnInformation(SDDS_DATASET *SDDS_dataset, char *field_name, void *memory, int32_t mode,...)
Retrieves information about a specified column in the SDDS dataset.
Definition SDDS_info.c:41
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.
int32_t SDDS_NumberOfErrors()
Retrieves the number of errors recorded by SDDS library routines.
Definition SDDS_utils.c:304
int32_t SDDS_StringIsBlank(char *s)
Checks if a string is blank (contains only whitespace characters).
int32_t SDDS_CopyString(char **target, const char *source)
Copies a source string to a target string with memory allocation.
Definition SDDS_utils.c:856
#define SDDS_STRING
Identifier for the string data type.
Definition SDDStypes.h:85
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37

◆ create_fft_frequency_column()

long create_fft_frequency_column ( SDDS_DATASET * SDDSout,
SDDS_DATASET * SDDSin,
char * timeName,
char * freqUnits,
long inverse )

Definition at line 875 of file sddsfft.c.

875 {
876 char s[SDDS_MAXLINE];
877 char *timeSymbol;
878 char *description;
879
880 if (SDDS_GetColumnInformation(SDDSin, "symbol", &timeSymbol, SDDS_GET_BY_NAME, timeName) != SDDS_STRING)
881 return 0;
882 if (!timeSymbol || SDDS_StringIsBlank(timeSymbol))
883 SDDS_CopyString(&timeSymbol, timeName);
884
885 sprintf(s, "Frequency for %s", timeSymbol);
886 SDDS_CopyString(&description, s);
887 if (!inverse) {
888 if (SDDS_DefineColumn(SDDSout, "f", NULL, freqUnits, description, NULL, SDDS_DOUBLE, 0) < 0) {
889 free(timeSymbol);
890 free(description);
891 return 0;
892 }
893 } else {
894 sprintf(s, "inverse for %s", timeSymbol);
895 SDDS_CopyString(&description, s);
896 if (SDDS_DefineColumn(SDDSout, "t", NULL, freqUnits, description, NULL, SDDS_DOUBLE, 0) < 0) {
897 free(timeSymbol);
898 free(description);
899 return 0;
900 }
901 }
902 free(timeSymbol);
903 free(description);
904 return 1;
905}

◆ expandComplexColumnPairNames()

long expandComplexColumnPairNames ( SDDS_DATASET * SDDSin,
char ** name,
char *** realName,
char *** imagName,
long names,
char ** excludeName,
long excludeNames,
long typeMode,
long typeValue )

Definition at line 1143 of file sddsfft.c.

1143 {
1144 long i, j, k, realNames, imagNames, names2;
1145 char **realName1, **imagName1, **realName2, **imagName2;
1146 char *realPattern, *imagPattern = NULL;
1147 long longest=0;
1148
1149 if (!names || !name)
1150 return 0;
1151 realName1 = imagName1 = realName2 = imagName2 = NULL;
1152 realNames = imagNames = names2 = 0;
1153 for (i = 0; i < names; i++) {
1154 if (strlen(name[i]) > longest)
1155 longest = strlen(name[i]);
1156 }
1157 longest += 10;
1158 if (!(realPattern = SDDS_Malloc(sizeof(*realPattern) * longest)) || !(imagPattern = SDDS_Malloc(sizeof(*imagPattern) * longest)))
1159 SDDS_Bomb("memory allocation failure");
1160
1161 for (i = 0; i < names; i++) {
1162 for (j = 0; j < 2; j++) {
1163 if (j == 0) {
1164 sprintf(realPattern, "Real%s", name[i]);
1165 sprintf(imagPattern, "Imag%s", name[i]);
1166 } else {
1167 sprintf(realPattern, "%sReal", name[i]);
1168 sprintf(imagPattern, "%sImag", name[i]);
1169 }
1170 switch (typeMode) {
1171 case FIND_ANY_TYPE:
1172 case FIND_NUMERIC_TYPE:
1173 case FIND_INTEGER_TYPE:
1174 case FIND_FLOATING_TYPE:
1175 realNames = SDDS_MatchColumns(SDDSin, &realName1, SDDS_MATCH_STRING, typeMode, realPattern, SDDS_0_PREVIOUS | SDDS_OR);
1176 imagNames = SDDS_MatchColumns(SDDSin, &imagName1, SDDS_MATCH_STRING, typeMode, imagPattern, SDDS_0_PREVIOUS | SDDS_OR);
1177 break;
1178 case FIND_SPECIFIED_TYPE:
1179 if (!SDDS_VALID_TYPE(typeValue))
1180 SDDS_Bomb("invalid type value in expandColumnPairNames");
1181 realNames = SDDS_MatchColumns(SDDSin, &realName1, SDDS_MATCH_STRING, typeMode, typeValue, realPattern, SDDS_0_PREVIOUS | SDDS_OR);
1182 imagNames = SDDS_MatchColumns(SDDSin, &imagName1, SDDS_MATCH_STRING, typeMode, typeValue, imagPattern, SDDS_0_PREVIOUS | SDDS_OR);
1183 break;
1184 default:
1185 SDDS_Bomb("invalid typeMode in expandColumnPairNames");
1186 exit(EXIT_FAILURE);
1187 break;
1188 }
1189 if (realNames == 0)
1190 continue;
1191 if (realNames == -1 || imagNames == -1) {
1192 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
1193 SDDS_Bomb("unable to perform column name match in expandColumnPairNames");
1194 }
1195 if (realNames != imagNames)
1196 SDDS_Bomb("found different number of real and imaginary columns");
1197 if (excludeNames) {
1198 for (j = 0; j < excludeNames; j++)
1199 for (k = 0; k < realNames; k++)
1200 if (wild_match(realName1[k], excludeName[j])) {
1201 free(realName1[k]);
1202 free(imagName1[k]);
1203 imagName1[k] = realName1[k] = NULL;
1204 }
1205 }
1206 moveToStringArrayComplex(&realName2, &imagName2, &names2, realName1, imagName1, realNames);
1207 free(realName1);
1208 free(imagName1);
1209 }
1210 }
1211 free(realPattern);
1212 free(imagPattern);
1213 if (names2 == 0)
1214 return 0;
1215 *realName = realName2;
1216 *imagName = imagName2;
1217 return names2;
1218}
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:432
void * SDDS_Malloc(size_t size)
Allocates memory of a specified size.
Definition SDDS_utils.c:639
int32_t SDDS_MatchColumns(SDDS_DATASET *SDDS_dataset, char ***nameReturn, int32_t matchMode, int32_t typeMode,...)
Matches and retrieves column names from an SDDS dataset based on specified criteria.
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:342
#define SDDS_VALID_TYPE(type)
Validates whether the given type identifier is within the defined range of SDDS types.
Definition SDDStypes.h:149
int wild_match(char *string, char *template)
Determine whether one string is a wildcard match for another.
Definition wild_match.c:49

◆ greatestProductOfSmallPrimes()

int64_t greatestProductOfSmallPrimes ( int64_t rows)

Definition at line 311 of file SDDSutils.c.

313{
314 int64_t bestResult = 0, result, nPrimes;
315 static int64_t prime[MAXPRIMES] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
316 71, 73, 79, 83, 89, 97};
317
318 for (nPrimes = 1; nPrimes <= MAXPRIMES; nPrimes++) {
319 if ((result = greatestProductOfSmallPrimes1(rows, prime, nPrimes)) > bestResult &&
320 result <= rows)
321 bestResult = result;
322 }
323 if (bestResult == 0)
324 SDDS_Bomb("couldn't find acceptable number of rows for truncation/padding");
325 return bestResult;
326}

◆ main()

int main ( int argc,
char ** argv )

Definition at line 203 of file sddsfft.c.

203 {
204 int iArg;
205 char *freqUnits;
206 char *indepQuantity, **depenQuantity, **exclude, **realQuan = NULL, **imagQuan = NULL;
207 long depenQuantities, excludes;
208 char *input, *output;
209 long i, readCode, noWarnings, complexInput, inverse, spectrumFoldParExist = 0;
210 int64_t rows, rowsToUse, sampleInterval;
211 int32_t spectrumFolded = 0, page = 0;
212 unsigned long flags, pipeFlags, complexInputFlags = 0, fullOutputFlags = 0, majorOrderFlag;
213 long windowType = -1;
214 SCANNED_ARG *scanned;
215 SDDS_DATASET SDDSin, SDDSout;
216 double *tdata, rintegCutOffFreq, unwrapLimit = 0;
217 long padFactor, correctWindowEffects = 0;
218 short columnMajorOrder = -1;
219
221 argc = scanargs(&scanned, argc, argv);
222 if (argc < 3 || argc > (3 + N_OPTIONS)) {
223 fprintf(stderr, "%s%s", USAGE1, USAGE2);
224 exit(EXIT_FAILURE);
225 /* bomb(NULL, USAGE); */
226 }
227 rintegCutOffFreq = 0;
228 output = input = NULL;
229 flags = pipeFlags = excludes = complexInput = inverse = 0;
230 sampleInterval = 1;
231 indepQuantity = NULL;
232 depenQuantity = exclude = NULL;
233 depenQuantities = 0;
234 noWarnings = 0;
235 padFactor = 0;
236 for (iArg = 1; iArg < argc; iArg++) {
237 if (scanned[iArg].arg_type == OPTION) {
238 /* process options here */
239 switch (match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
240 case SET_NORMALIZE:
241 flags |= FL_NORMALIZE;
242 break;
243 case SET_WINDOW:
244 if (scanned[iArg].n_items != 1) {
245 if ((i = match_string(scanned[iArg].list[1], window_type, N_WINDOW_TYPES, 0)) < 0)
246 SDDS_Bomb("unknown window type");
247 windowType = i;
248 if (scanned[iArg].n_items > 2) {
249 if (strncmp(scanned[iArg].list[2], "correct", strlen(scanned[iArg].list[2])) == 0)
250 correctWindowEffects = 1;
251 else
252 SDDS_Bomb("invalid -window syntax");
253 }
254 } else
255 windowType = 0;
256 break;
257 case SET_PADWITHZEROES:
258 flags |= FL_PADWITHZEROES;
259 if (scanned[iArg].n_items != 1) {
260 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%ld", &padFactor) != 1 || padFactor < 1)
261 SDDS_Bomb("invalid -padwithzeroes syntax");
262 }
263 break;
264 case SET_TRUNCATE:
265 flags |= FL_TRUNCATE;
266 break;
267 case SET_SUPPRESSAVERAGE:
268 flags |= FL_SUPPRESSAVERAGE;
269 break;
270 case SET_SAMPLEINTERVAL:
271 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%" SCNd64, &sampleInterval) != 1 || sampleInterval <= 0)
272 SDDS_Bomb("invalid -sampleinterval syntax");
273 break;
274 case SET_COLUMNS:
275 if (indepQuantity)
276 SDDS_Bomb("only one -columns option may be given");
277 if (scanned[iArg].n_items < 2)
278 SDDS_Bomb("invalid -columns syntax");
279 indepQuantity = scanned[iArg].list[1];
280 if (scanned[iArg].n_items >= 2) {
281 depenQuantity = tmalloc(sizeof(*depenQuantity) * (depenQuantities = scanned[iArg].n_items - 2));
282 for (i = 0; i < depenQuantities; i++)
283 depenQuantity[i] = scanned[iArg].list[i + 2];
284 }
285 break;
286 case SET_FULLOUTPUT:
287 flags |= FL_FULLOUTPUT;
288 if (scanned[iArg].n_items >= 2) {
289 scanned[iArg].n_items--;
290 if (!scanItemList(&fullOutputFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "folded", -1, NULL, 0, FL_FULLOUTPUT_FOLDED, "unfolded", -1, NULL, 0, FL_FULLOUTPUT_UNFOLDED, "unwrapLimit", SDDS_DOUBLE, &unwrapLimit, 0, FL_UNWRAP_PHASE, NULL))
291 SDDS_Bomb("Invalid -fullOutput syntax");
292 scanned[iArg].n_items++;
293 if (fullOutputFlags & FL_FULLOUTPUT_UNFOLDED)
294 flags |= FL_FULLOUTPUT_UNFOLDED;
295 else
296 flags |= FL_FULLOUTPUT_FOLDED;
297 if (fullOutputFlags & FL_UNWRAP_PHASE)
298 flags |= FL_UNWRAP_PHASE;
299 }
300 break;
301 case SET_PIPE:
302 if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
303 SDDS_Bomb("invalid -pipe syntax");
304 break;
305 case SET_PSDOUTPUT:
306 if (scanned[iArg].n_items -= 1) {
307 unsigned long tmpFlags;
308 if (strchr(scanned[iArg].list[1], '=') <= 0) {
309 if (!scanItemList(&tmpFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "integrated", -1, NULL, 0, FL_PSDINTEGOUTPUT, "rintegrated", -1, NULL, 0, FL_PSDRINTEGOUTPUT, "plain", -1, NULL, 0, FL_PSDOUTPUT, NULL))
310 SDDS_Bomb("invalid -psdOutput syntax");
311 } else {
312 if (!scanItemList(&tmpFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "integrated", -1, NULL, 0, FL_PSDINTEGOUTPUT, "rintegrated", SDDS_DOUBLE, &rintegCutOffFreq, 0, FL_PSDRINTEGOUTPUT, "plain", -1, NULL, 0, FL_PSDOUTPUT, NULL))
313 SDDS_Bomb("invalid -psdOutput syntax");
314 }
315 flags |= tmpFlags;
316 } else
317 flags |= FL_PSDOUTPUT;
318 if (flags & FL_PSDINTEGOUTPUT && flags & FL_PSDRINTEGOUTPUT)
319 SDDS_Bomb("invalid -psdOutput syntax: give only one of integrated or rintegrated");
320 break;
321 case SET_EXCLUDE:
322 if (scanned[iArg].n_items < 2)
323 SDDS_Bomb("invalid -exclude syntax");
324 moveToStringArray(&exclude, &excludes, scanned[iArg].list + 1, scanned[iArg].n_items - 1);
325 break;
326 case SET_NOWARNINGS:
327 noWarnings = 1;
328 break;
329 case SET_COMPLEXINPUT:
330 complexInput = 1;
331 if (scanned[iArg].n_items == 2) {
332 scanned[iArg].n_items--;
333 if (!scanItemList(&complexInputFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "folded", -1, NULL, 0, FL_COMPLEXINPUT_FOLDED, "unfolded", -1, NULL, 0, FL_COMPLEXINPUT_UNFOLDED, NULL))
334 SDDS_Bomb("Invalid -complexInput syntax");
335 scanned[iArg].n_items++;
336 }
337 break;
338 case SET_INVERSE:
339 inverse = 1;
340 break;
341 case SET_MAJOR_ORDER:
342 majorOrderFlag = 0;
343 scanned[iArg].n_items--;
344 if (scanned[iArg].n_items > 0 && (!scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER, "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
345 SDDS_Bomb("invalid -majorOrder syntax/values");
346 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
347 columnMajorOrder = 1;
348 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
349 columnMajorOrder = 0;
350 break;
351 default:
352 fprintf(stderr, "error: unknown/ambiguous option: %s\n", scanned[iArg].list[0]);
353 exit(EXIT_FAILURE);
354 break;
355 }
356 } else {
357 if (!input)
358 input = scanned[iArg].list[0];
359 else if (!output)
360 output = scanned[iArg].list[0];
361 else
362 SDDS_Bomb("too many filenames seen");
363 }
364 }
365 if (!complexInput) {
366 if (!noWarnings && inverse)
367 fprintf(stderr, "Warning: the inverse option is ignored since it only works with -complexInput.\n");
368 inverse = 0;
369 }
370 if (!noWarnings && inverse && flags & FL_FULLOUTPUT_FOLDED)
371 fprintf(stderr, "Warning: the -inverse -fullOutput=folded will be changed to -inverse -fullOutput=unfolded.\n");
372
373 processFilenames("sddsfft", &input, &output, pipeFlags, 0, NULL);
374
375 if (!indepQuantity)
376 SDDS_Bomb("Supply the independent quantity name with the -columns option.");
377
378 if (flags & FL_TRUNCATE && flags & FL_PADWITHZEROES)
379 SDDS_Bomb("Specify only one of -padwithzeroes and -truncate.");
380
381 if (!SDDS_InitializeInput(&SDDSin, input))
382 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
383
384 if (SDDS_CheckColumn(&SDDSin, indepQuantity, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OKAY)
385 exit(EXIT_FAILURE);
386
387 excludes = appendToStringArray(&exclude, excludes, indepQuantity);
388 if (!depenQuantities)
389 depenQuantities = appendToStringArray(&depenQuantity, depenQuantities, "*");
390
391 if (!complexInput) {
392 if ((depenQuantities = expandColumnPairNames(&SDDSin, &depenQuantity, NULL, depenQuantities, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
393 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
394 SDDS_Bomb("No quantities selected to FFT.");
395 }
396 } else {
397 if ((depenQuantities = expandComplexColumnPairNames(&SDDSin, depenQuantity, &realQuan, &imagQuan, depenQuantities, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
398 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
399 SDDS_Bomb("No quantities selected to FFT.");
400 }
401 }
402
403#if 0
404 fprintf(stderr, "%ld dependent quantities:\n", depenQuantities);
405 for (i = 0; i < depenQuantities; i++)
406 fprintf(stderr, " %s\n", depenQuantity[i]);
407#endif
408
409 if (!(freqUnits = makeFrequencyUnits(&SDDSin, indepQuantity)) ||
410 !SDDS_InitializeOutput(&SDDSout, SDDS_BINARY, 0, NULL, "sddsfft output", output) ||
411 !create_fft_frequency_column(&SDDSout, &SDDSin, indepQuantity, freqUnits, inverse) ||
412 SDDS_DefineParameter(&SDDSout, "fftFrequencies", NULL, NULL, NULL, NULL, SDDS_LONG, NULL) < 0 ||
413 SDDS_DefineParameter(&SDDSout, "fftFrequencySpacing", "$gD$rf", freqUnits, NULL, NULL, SDDS_DOUBLE, NULL) < 0)
414 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
415 if (columnMajorOrder != -1)
416 SDDSout.layout.data_mode.column_major = columnMajorOrder;
417 else
418 SDDSout.layout.data_mode.column_major = SDDSin.layout.data_mode.column_major;
419
420 if (flags & FL_FULLOUTPUT && SDDS_DefineParameter(&SDDSout, "SpectrumFolded", NULL, NULL, NULL, NULL, SDDS_LONG, NULL) < 0)
421 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
422 if (complexInput) {
423 if (!complexInputFlags) {
424 if (SDDS_CheckParameter(&SDDSin, "SpectrumFolded", NULL, SDDS_LONG, NULL) == SDDS_CHECK_OK)
425 spectrumFoldParExist = 1;
426 } else if (complexInputFlags & FL_COMPLEXINPUT_UNFOLDED)
427 flags |= FL_COMPLEXINPUT_UNFOLDED;
428 else
429 flags |= FL_COMPLEXINPUT_FOLDED;
430 }
431 for (i = 0; i < depenQuantities; i++) {
432 if (!complexInput)
433 create_fft_columns(&SDDSout, &SDDSin, depenQuantity[i], indepQuantity, freqUnits,
434 flags & FL_FULLOUTPUT,
435 flags & (FL_PSDOUTPUT + FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT),
436 0, inverse, flags & FL_UNWRAP_PHASE);
437 else
438 create_fft_columns(&SDDSout, &SDDSin, realQuan[i], indepQuantity, freqUnits,
439 flags & FL_FULLOUTPUT,
440 flags & (FL_PSDOUTPUT + FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT),
441 1, inverse, flags & FL_UNWRAP_PHASE);
442 }
443
444 if (!SDDS_TransferAllParameterDefinitions(&SDDSout, &SDDSin, SDDS_TRANSFER_KEEPOLD) ||
445 !SDDS_WriteLayout(&SDDSout))
446 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
447
448 while ((readCode = SDDS_ReadPage(&SDDSin)) > 0) {
449 page++;
450 if ((rows = SDDS_CountRowsOfInterest(&SDDSin)) < 0)
451 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
452 if (page == 1 && spectrumFoldParExist) {
453 if (!SDDS_GetParameterAsLong(&SDDSin, "SpectrumFolded", &spectrumFolded))
454 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
455 if (spectrumFolded)
456 flags |= FL_COMPLEXINPUT_FOLDED;
457 else
458 flags |= FL_COMPLEXINPUT_UNFOLDED;
459 }
460 if (rows) {
461 int64_t primeRows, pow2Rows;
462 rowsToUse = rows;
463 primeRows = greatestProductOfSmallPrimes(rows);
464 if (rows != primeRows || padFactor) {
465 if (flags & FL_PADWITHZEROES) {
466 pow2Rows = ipow(2., ((int64_t)(log((double)rows) / log(2.0F))) + (padFactor ? padFactor : 1.0));
467 if ((primeRows = greatestProductOfSmallPrimes(pow2Rows)) > rows)
468 rowsToUse = primeRows;
469 else
470 rowsToUse = pow2Rows;
471 } else if (flags & FL_TRUNCATE)
472 rowsToUse = greatestProductOfSmallPrimes(rows);
473 else if (largest_prime_factor(rows) > 1000 && !noWarnings)
474 fputs("Warning: number of points has large prime factors.\nThis could take a very long time.\nConsider using the -truncate option.\n", stderr);
475 }
476 if (!SDDS_StartPage(&SDDSout, 2 * rowsToUse + 2) || !SDDS_CopyParameters(&SDDSout, &SDDSin))
477 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
478 if (!(tdata = SDDS_GetColumnInDoubles(&SDDSin, indepQuantity)))
479 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
480 for (i = 0; i < depenQuantities; i++)
481 if (!process_data(&SDDSout, &SDDSin, tdata, rows, rowsToUse,
482 complexInput ? realQuan[i] : depenQuantity[i],
483 complexInput ? imagQuan[i] : NULL,
484 flags | (i == 0 ? FL_MAKEFREQDATA : 0),
485 windowType, sampleInterval, correctWindowEffects, inverse,
486 rintegCutOffFreq, unwrapLimit)) {
487 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
488 exit(EXIT_FAILURE);
489 }
490 free(tdata);
491 } else {
492 if (!SDDS_StartPage(&SDDSout, 0) || !SDDS_CopyParameters(&SDDSout, &SDDSin))
493 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
494 }
495 if (!SDDS_WritePage(&SDDSout))
496 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
497 }
498
499 if (!SDDS_Terminate(&SDDSin)) {
500 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
501 exit(EXIT_FAILURE);
502 }
503 if (!SDDS_Terminate(&SDDSout)) {
504 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
505 exit(EXIT_FAILURE);
506 }
507
508 return EXIT_SUCCESS;
509}
int32_t SDDS_CopyParameters(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
Definition SDDS_copy.c:286
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
int32_t * SDDS_GetParameterAsLong(SDDS_DATASET *SDDS_dataset, char *parameter_name, int32_t *memory)
Retrieves the value of a specified parameter as a 32-bit integer from the current data table of a dat...
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)
int32_t SDDS_InitializeOutput(SDDS_DATASET *SDDS_dataset, int32_t data_mode, int32_t lines_per_row, const char *description, const char *contents, const char *filename)
Initializes the SDDS output dataset.
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_DefineParameter(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, char *fixed_value)
Defines a data parameter with a fixed string value.
int32_t SDDS_TransferAllParameterDefinitions(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source, uint32_t mode)
Transfers all parameter definitions from a source dataset to a target dataset.
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
int32_t SDDS_CheckParameter(SDDS_DATASET *SDDS_dataset, char *name, char *units, int32_t type, FILE *fp_message)
Checks if a parameter exists in the SDDS dataset with the specified name, units, and type.
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
#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
int64_t largest_prime_factor(int64_t number)
Find the largest prime factor of a number.
Definition factorize.c:83
double ipow(const double x, const int64_t p)
Compute x raised to the power p (x^p).
Definition ipow.c:33
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
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.

◆ makeFrequencyUnits()

char * makeFrequencyUnits ( SDDS_DATASET * SDDSin,
char * indepName )

Definition at line 269 of file SDDSutils.c.

269 {
270 char *timeUnits;
271 char *units;
272 long reciprocal = 0, end;
273
274 if (SDDS_GetColumnInformation(SDDSin, "units", &timeUnits, SDDS_GET_BY_NAME, indepName) != SDDS_STRING)
275 return 0;
276 if (timeUnits) {
277 while (1) {
278 end = strlen(timeUnits) - 1;
279 if (timeUnits[0] == '(' && timeUnits[end] == ')') {
280 timeUnits[end] = 0;
281 strslide(timeUnits, 1);
282 } else if (timeUnits[0] == '1' && timeUnits[1] == '/' && timeUnits[2] == '(' && timeUnits[end] == ')') {
283 timeUnits[end] = 0;
284 strslide(timeUnits, 3);
285 reciprocal = !reciprocal;
286 } else
287 break;
288 }
289 }
290 if (!timeUnits || SDDS_StringIsBlank(timeUnits)) {
291 units = tmalloc(sizeof(*units) * 1);
292 units[0] = 0;
293 return units;
294 }
295
296 if (reciprocal) {
297 if (!SDDS_CopyString(&units, timeUnits))
298 return NULL;
299 return units;
300 }
301
302 units = tmalloc(sizeof(*units) * (strlen(timeUnits) + 5));
303 if (strchr(timeUnits, ' '))
304 sprintf(units, "1/(%s)", timeUnits);
305 else
306 sprintf(units, "1/%s", timeUnits);
307 return units;
308}
char * strslide(char *s, long distance)
Slides character data within a string by a specified distance.
Definition strslide.c:32

◆ moveToStringArrayComplex()

void moveToStringArrayComplex ( char *** targetReal,
char *** targetImag,
long * targets,
char ** sourceReal,
char ** sourceImag,
long sources )

Definition at line 1220 of file sddsfft.c.

1220 {
1221 long i, j;
1222 if (!sources)
1223 return;
1224 if (!(*targetReal = SDDS_Realloc(*targetReal, sizeof(**targetReal) * (*targets + sources))) ||
1225 !(*targetImag = SDDS_Realloc(*targetImag, sizeof(**targetImag) * (*targets + sources))))
1226 SDDS_Bomb("memory allocation failure");
1227 for (i = 0; i < sources; i++) {
1228 if (sourceReal[i] == NULL || sourceImag[i] == NULL)
1229 continue;
1230 for (j = 0; j < *targets; j++)
1231 if (strcmp(sourceReal[i], (*targetReal)[j]) == 0)
1232 break;
1233 if (j == *targets) {
1234 (*targetReal)[j] = sourceReal[i];
1235 (*targetImag)[j] = sourceImag[i];
1236 *targets += 1;
1237 }
1238 }
1239}
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:677

◆ process_data()

long process_data ( SDDS_DATASET * SDDSout,
SDDS_DATASET * SDDSin,
double * tdata,
int64_t rows,
int64_t rowsToUse,
char * depenQuantity,
char * depenQuantity2,
unsigned long flags,
long windowType,
int64_t sampleInterval,
long correctWindowEffects,
long inverse,
double rintegCutOffFreq,
double unwrapLimit )

Definition at line 513 of file sddsfft.c.

515 {
516 long offset, index, unfold = 0;
517 int64_t n_freq, i, fftrows = 0;
518 double r, r1, r2, length, factor, df, min, max, delta;
519 double *real, *imag, *magData, *arg = NULL, *real_imag, *data, *psd = NULL, *psdInteg = NULL, *psdIntegPower = NULL, *unwrapArg = NULL, phase_correction = 0;
520 double dtf_real, dtf_imag, t0;
521 char s[256];
522 double *fdata, *imagData = NULL;
523 double *tDataStore = NULL;
524 double windowCorrectionFactor = 0;
525
526 if (!(data = SDDS_GetColumnInDoubles(SDDSin, depenQuantity)))
527 return 0;
528 if (imagQuantity && !(imagData = SDDS_GetColumnInDoubles(SDDSin, imagQuantity)))
529 return 0;
530 if (flags & FL_SUPPRESSAVERAGE) {
531 compute_average(&r, data, rows);
532 for (i = 0; i < rows; i++)
533 data[i] -= r;
534 if (imagData) {
535 compute_average(&r, imagData, rows);
536 for (i = 0; i < rows; i++)
537 imagData[i] -= r;
538 }
539 }
540 if (rows < rowsToUse) {
541 /* pad with zeroes */
542 tDataStore = tmalloc(sizeof(*tDataStore) * rowsToUse);
543 memcpy((char *)tDataStore, (char *)tdata, rows * sizeof(*tdata));
544 if (!(data = SDDS_Realloc(data, sizeof(*data) * rowsToUse)))
545 SDDS_Bomb("memory allocation failure");
546 if (imagData && !(imagData = SDDS_Realloc(imagData, sizeof(*imagData) * rowsToUse)))
547 length = ((double)rows) * (tdata[rows - 1] - tdata[0]) / ((double)rows - 1.0);
548 else
549 length = tdata[rows - 1] - tdata[0];
550 for (i = rows; i < rowsToUse; i++) {
551 tDataStore[i] = tDataStore[i - 1] + length / ((double)rows - 1);
552 data[i] = 0;
553 }
554 if (imagData)
555 for (i = rows; i < rowsToUse; i++)
556 imagData[i] = 0;
557 tdata = tDataStore;
558 }
559 rows = rowsToUse; /* results in truncation if rows>rowsToUse */
560 windowCorrectionFactor = 0;
561 switch (windowType) {
562 case WINDOW_HANNING:
563 r = PIx2 / (rows - 1);
564 for (i = 0; i < rows; i++) {
565 factor = (1 - cos(i * r)) / 2;
566 data[i] *= factor;
567 windowCorrectionFactor += sqr(factor);
568 if (imagData)
569 imagData[i] *= factor;
570 }
571 break;
572 case WINDOW_HAMMING:
573 r = PIx2 / (rows - 1);
574 for (i = 0; i < rows; i++) {
575 factor = 0.54 - 0.46 * cos(i * r);
576 data[i] *= factor;
577 windowCorrectionFactor += sqr(factor);
578 if (imagData)
579 imagData[i] *= factor;
580 }
581 if (imagData)
582 for (i = 0; i < rows; i++)
583 imagData[i] *= (1 - cos(i * r)) / 2;
584 break;
585 case WINDOW_WELCH:
586 r1 = (rows - 1) / 2.0;
587 r2 = sqr((rows + 1) / 2.0);
588 for (i = 0; i < rows; i++) {
589 factor = 1 - sqr(i - r1) / r2;
590 data[i] *= factor;
591 windowCorrectionFactor += sqr(factor);
592 if (imagData)
593 imagData[i] *= factor;
594 }
595 break;
596 case WINDOW_PARZEN:
597 r = (rows - 1) / 2.0;
598 for (i = 0; i < rows; i++) {
599 factor = 1 - FABS((i - r) / r);
600 data[i] *= factor;
601 windowCorrectionFactor += sqr(factor);
602 if (imagData)
603 imagData[i] *= factor;
604 }
605 break;
606 case WINDOW_FLATTOP:
607 for (i = 0; i < rows; i++) {
608 r = i * PIx2 / (rows - 1);
609 factor = 1 - 1.93 * cos(r) + 1.29 * cos(2 * r) - 0.388 * cos(3 * r) + 0.032 * cos(4 * r);
610 data[i] *= factor;
611 windowCorrectionFactor += sqr(factor);
612 if (imagData)
613 imagData[i] *= factor;
614 }
615 break;
616 case WINDOW_GAUSSIAN:
617 for (i = 0; i < rows; i++) {
618 r = sqr((i - (rows - 1) / 2.) / (0.4 * (rows - 1) / 2.)) / 2;
619 factor = exp(-r);
620 data[i] *= factor;
621 windowCorrectionFactor += sqr(factor);
622 if (imagData)
623 imagData[i] *= factor;
624 }
625 break;
626 case WINDOW_NONE:
627 default:
628 windowCorrectionFactor = 1;
629 break;
630 }
631
632 if (correctWindowEffects) {
633 /* Add correction factor to make the integrated PSD come out right. */
634 windowCorrectionFactor = 1 / sqrt(windowCorrectionFactor / rows);
635 for (i = 0; i < rows; i++)
636 data[i] *= windowCorrectionFactor;
637 if (imagData)
638 imagData[i] *= windowCorrectionFactor;
639 }
640 if (imagData && flags & FL_COMPLEXINPUT_FOLDED) {
641 double min, max, max1;
642 data = SDDS_Realloc(data, sizeof(*data) * rows * 2);
643 imagData = SDDS_Realloc(imagData, sizeof(*data) * rows * 2);
644 find_min_max(&min, &max, data, rows);
645 if (fabs(min) > fabs(max))
646 max1 = fabs(min);
647 else
648 max1 = fabs(max);
649 find_min_max(&min, &max, imagData, rows);
650 if (fabs(min) > max1)
651 max1 = fabs(min);
652 if (fabs(max) > max1)
653 max1 = fabs(max);
654 if (fabs(imagData[rows - 1]) / max1 < 1.0e-15) {
655 fftrows = 2 * (rows - 1);
656 for (i = 1; i < rows - 1; i++) {
657 data[i] = data[i] / 2.0;
658 imagData[i] = imagData[i] / 2.0;
659 }
660 for (i = 1; i < rows - 1; i++) {
661 data[rows - 1 + i] = data[rows - 1 - i];
662 imagData[rows - 1 + i] = -imagData[rows - 1 - i];
663 }
664 length = (tdata[rows - 1] - tdata[0]) * 2.0;
665 } else {
666 fftrows = 2 * (rows - 1) + 1;
667 for (i = 1; i < rows; i++) {
668 data[i] = data[i] / 2.0;
669 imagData[i] = imagData[i] / 2.0;
670 }
671 for (i = 0; i < rows - 1; i++) {
672 data[rows + i] = data[rows - 1 - i];
673 imagData[rows + i] = -imagData[rows - 1 - i];
674 }
675 length = ((double)fftrows) * (tdata[rows - 1] - tdata[0]) / ((double)fftrows - 1.0) * 2;
676 }
677 } else {
678 fftrows = rows;
679 length = ((double)rows) * (tdata[rows - 1] - tdata[0]) / ((double)rows - 1.0);
680 }
681 /* compute FFT */
682
683 real_imag = tmalloc(sizeof(double) * (2 * fftrows + 2));
684 for (i = 0; i < fftrows; i++) {
685 real_imag[2 * i] = data[i];
686 if (imagData)
687 real_imag[2 * i + 1] = imagData[i];
688 else
689 real_imag[2 * i + 1] = 0;
690 }
691 if (!inverse) {
692 complexFFT(real_imag, fftrows, 0);
693 if (flags & FL_FULLOUTPUT_UNFOLDED) {
694 n_freq = fftrows;
695 unfold = 1;
696 } else if (flags & FL_FULLOUTPUT_FOLDED)
697 n_freq = fftrows / 2 + 1;
698 else if (!imagData)
699 n_freq = fftrows / 2 + 1;
700 else
701 n_freq = fftrows + 1;
702 } else {
703 complexFFT(real_imag, fftrows, INVERSE_FFT);
704 n_freq = fftrows;
705 }
706 /* calculate factor for converting k to f or omega */
707 /* length is assumed the length of period, not the total length of the data file */
708
709 /* length = ((double)rows)*(tdata[rows-1]-tdata[0])/((double)rows-1.0); */
710 t0 = tdata[0];
711 df = factor = 1.0 / length;
712
713 /* convert into amplitudes and frequencies, adding phase factor for t[0]!=0 */
714 real = tmalloc(sizeof(double) * n_freq);
715 imag = tmalloc(sizeof(double) * n_freq);
716 fdata = tmalloc(sizeof(double) * n_freq);
717 magData = tmalloc(sizeof(double) * n_freq);
718 if (flags & FL_PSDOUTPUT || flags & (FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT)) {
719 psd = tmalloc(sizeof(*psd) * n_freq);
720 if (flags & (FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT)) {
721 psdInteg = tmalloc(sizeof(*psdInteg) * n_freq);
722 psdIntegPower = tmalloc(sizeof(*psdIntegPower) * n_freq);
723 }
724 }
725
726 for (i = 0; i < n_freq; i++) {
727 fdata[i] = i * df;
728 dtf_real = cos(-2 * PI * fdata[i] * t0);
729 dtf_imag = sin(-2 * PI * fdata[i] * t0);
730 if (psd)
731 psd[i] = (sqr(real_imag[2 * i]) + sqr(real_imag[2 * i + 1])) / df;
732 if (!imagData && i != 0 && !(i == (n_freq - 1) && rows % 2 == 0)) {
733 /* This is not the DC or Nyquist term, so
734 multiply by 2 to account for amplitude in
735 negative frequencies
736 */
737 if (!unfold) {
738 real_imag[2 * i] *= 2;
739 real_imag[2 * i + 1] *= 2;
740 }
741 if (psd)
742 psd[i] *= 2; /* 2 really is what I want--not 4 */
743 }
744 real[i] = real_imag[2 * i] * dtf_real - real_imag[2 * i + 1] * dtf_imag;
745 imag[i] = real_imag[2 * i + 1] * dtf_real + real_imag[2 * i] * dtf_imag;
746 magData[i] = sqrt(sqr(real[i]) + sqr(imag[i]));
747 }
748
749 if (psdInteg) {
750 if (flags & FL_PSDINTEGOUTPUT) {
751 psdIntegPower[0] = 0;
752 for (i = 1; i < n_freq; i++)
753 psdIntegPower[i] = psdIntegPower[i - 1] + (psd[i - 1] + psd[i]) * df / 2.;
754 for (i = 0; i < n_freq; i++)
755 psdInteg[i] = sqrt(psdIntegPower[i]);
756 } else {
757 psdIntegPower[n_freq - 1] = 0;
758 for (i = n_freq - 2; i >= 0; i--) {
759 if (rintegCutOffFreq == 0 || fdata[i] <= rintegCutOffFreq)
760 psdIntegPower[i] = psdIntegPower[i + 1] + (psd[i + 1] + psd[i]) * df / 2.;
761 }
762 for (i = 0; i < n_freq; i++)
763 psdInteg[i] = sqrt(psdIntegPower[i]);
764 }
765 }
766
767 if (flags & FL_FULLOUTPUT) {
768 arg = tmalloc(sizeof(*arg) * n_freq);
769 for (i = 0; i < n_freq; i++) {
770 if (real[i] || imag[i])
771 arg[i] = 180.0 / PI * atan2(imag[i], real[i]);
772 else
773 arg[i] = 0;
774 }
775 }
776 if (flags & FL_UNWRAP_PHASE) {
777 find_min_max(&min, &max, magData, n_freq);
778 unwrapArg = tmalloc(sizeof(*unwrapArg) * n_freq);
779 phase_correction = 0;
780 for (i = 0; i < n_freq; i++) {
781 if (i && magData[i] / max > unwrapLimit) {
782 delta = arg[i] - arg[i - 1];
783 if (delta < -180.0)
784 phase_correction += 360.0;
785 else if (delta > 180.0)
786 phase_correction -= 360.0;
787 }
788 unwrapArg[i] = arg[i] + phase_correction;
789 }
790 }
791
792 if (flags & FL_NORMALIZE) {
793 factor = -DBL_MAX;
794 for (i = 0; i < n_freq; i++)
795 if (magData[i] > factor)
796 factor = magData[i];
797 if (factor != -DBL_MAX)
798 for (i = 0; i < n_freq; i++) {
799 real[i] /= factor;
800 imag[i] /= factor;
801 magData[i] /= factor;
802 }
803 }
804 if (!inverse)
805 sprintf(s, "FFT%s", depenQuantity + (imagData ? 4 : 0));
806 else {
807 if (strncmp(depenQuantity, "FFT", 3) == 0)
808 sprintf(s, "%s", depenQuantity + 3);
809 else if (strncmp(depenQuantity, "RealFFT", 7) == 0)
810 sprintf(s, "%s", depenQuantity + 7);
811 else
812 sprintf(s, "%s", depenQuantity);
813 }
814
815 if ((index = SDDS_GetColumnIndex(SDDSout, s)) < 0)
816 return 0;
817
818 if (flags & FL_SUPPRESSAVERAGE) {
819 n_freq -= 1;
820 offset = 1;
821 } else
822 offset = 0;
823
824 if ((flags & FL_MAKEFREQDATA &&
825 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, fdata + offset, n_freq, 0)) ||
826 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, magData + offset, n_freq, index + fftOffset) ||
827 (flags & FL_FULLOUTPUT &&
828 (!SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, real + offset, n_freq, index + realOffset) ||
829 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, imag + offset, n_freq, index + imagOffset) ||
830 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, arg + offset, n_freq, index + argOffset))) ||
831 (flags & FL_PSDOUTPUT &&
832 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, psd + offset, n_freq, index + psdOffset)) ||
833 (flags & (FL_PSDINTEGOUTPUT + FL_PSDRINTEGOUTPUT) && (!SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, psdInteg + offset, n_freq, index + psdIntOffset) ||
834 !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, psdIntegPower + offset, n_freq, index + psdIntPowerOffset))) ||
835 (flags & FL_UNWRAP_PHASE && !SDDS_SetColumn(SDDSout, SDDS_SET_BY_INDEX, unwrapArg + offset, n_freq, index + unwrappedArgOffset)))
836 return 0;
837 if (sampleInterval > 0) {
838 int32_t *sample_row_flag;
839 sample_row_flag = calloc(sizeof(*sample_row_flag), n_freq);
840 for (i = 0; i < n_freq; i += sampleInterval)
841 sample_row_flag[i] = 1;
842 if (!SDDS_AssertRowFlags(SDDSout, SDDS_FLAG_ARRAY, sample_row_flag, n_freq)) {
843 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
844 return 0;
845 }
846 free(sample_row_flag);
847 }
848 if (!SDDS_SetParameters(SDDSout, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, "fftFrequencies", n_freq, "fftFrequencySpacing", df, NULL))
849 return 0;
850 if (flags & FL_FULLOUTPUT && !SDDS_SetParameters(SDDSout, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, "SpectrumFolded", flags & FL_FULLOUTPUT_UNFOLDED ? 0 : 1, NULL))
851 return 0;
852 free(data);
853 free(magData);
854 if (imagData)
855 free(imagData);
856 if (tDataStore)
857 free(tDataStore);
858 if (arg)
859 free(arg);
860 free(real_imag);
861 free(real);
862 free(imag);
863 free(fdata);
864 if (psd)
865 free(psd);
866 if (psdInteg)
867 free(psdInteg);
868 if (psdIntegPower)
869 free(psdIntegPower);
870 if (unwrapArg)
871 free(unwrapArg);
872 return 1;
873}
int32_t SDDS_SetParameters(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
int32_t SDDS_SetColumn(SDDS_DATASET *SDDS_dataset, int32_t mode, void *data, int64_t rows,...)
Sets the values for one data column in the current data table of an SDDS dataset.
int32_t SDDS_AssertRowFlags(SDDS_DATASET *SDDS_dataset, uint32_t mode,...)
Sets acceptance flags for rows based on specified criteria.
int32_t SDDS_GetColumnIndex(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the index of a named column in the SDDS dataset.
int find_min_max(double *min, double *max, double *list, int64_t n)
Finds the minimum and maximum values in a list of doubles.
Definition findMinMax.c:33
long compute_average(double *value, double *data, int64_t n)
Computes the average of an array of doubles.
Definition median.c:144

Variable Documentation

◆ argOffset

long argOffset
static

Definition at line 511 of file sddsfft.c.

◆ fftOffset

long fftOffset = -1
static

Definition at line 511 of file sddsfft.c.

◆ imagOffset

long imagOffset
static

Definition at line 511 of file sddsfft.c.

◆ option

char* option[N_OPTIONS]
Initial value:
= {
"window",
"normalize",
"padwithzeroes",
"truncate",
"suppressaverage",
"sampleinterval",
"columns",
"fulloutput",
"pipe",
"psdoutput",
"exclude",
"nowarnings",
"complexinput",
"inverse",
"majorOrder",
}

Definition at line 66 of file sddsfft.c.

66 {
67 "window",
68 "normalize",
69 "padwithzeroes",
70 "truncate",
71 "suppressaverage",
72 "sampleinterval",
73 "columns",
74 "fulloutput",
75 "pipe",
76 "psdoutput",
77 "exclude",
78 "nowarnings",
79 "complexinput",
80 "inverse",
81 "majorOrder",
82};

◆ psdIntOffset

long psdIntOffset
static

Definition at line 511 of file sddsfft.c.

◆ psdIntPowerOffset

long psdIntPowerOffset
static

Definition at line 511 of file sddsfft.c.

◆ psdOffset

long psdOffset
static

Definition at line 511 of file sddsfft.c.

◆ realOffset

long realOffset
static

Definition at line 511 of file sddsfft.c.

◆ unwrappedArgOffset

long unwrappedArgOffset = -1
static

Definition at line 511 of file sddsfft.c.

◆ USAGE1

char* USAGE1
static
Initial value:
=
"Usage:\n"
" sddsfft [<inputfile>] [<outputfile>]\n"
" [-pipe=[input][,output]]\n"
" [-columns=<indep-variable>[,<depen-quantity>[,...]]]\n"
" [-complexInput[=unfolded|folded]]\n"
" [-exclude=<depen-quantity>[,...]]\n"
" [-window[={hanning|welch|parzen|hamming|flattop|gaussian|none}[,correct]]]\n"
" [-sampleInterval=<number>]\n"
" [-normalize]\n"
" [-fullOutput[=unfolded|folded],unwrapLimit=<value>]\n"
" [-psdOutput[=plain][,{integrated|rintegrated[=<cutoff>]}]]\n"
" [-inverse]\n"
" [-padwithzeroes[=exponent] | -truncate]\n"
" [-suppressaverage]\n"
" [-noWarnings]\n"
" [-majorOrder=row|column]\n\n"

Definition at line 111 of file sddsfft.c.

◆ USAGE2

char* USAGE2
static

Definition at line 129 of file sddsfft.c.

◆ window_type

char* window_type[N_WINDOW_TYPES]
Initial value:
= {
"hanning", "welch", "parzen", "hamming", "flattop", "gaussian", "none"}

Definition at line 107 of file sddsfft.c.

107 {
108 "hanning", "welch", "parzen", "hamming", "flattop", "gaussian", "none"};