-
Notifications
You must be signed in to change notification settings - Fork 1
/
sacbutter.c
166 lines (137 loc) · 6.52 KB
/
sacbutter.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
/* Copyright (c) 2021 Nils Maercklin */
/************************************************************************
SACBUTTER - SAC Butterworth high-pass and/or low-pass filter
*************************************************************************
Author: Nils Maercklin,
RISSC, University of Naples, Italy, September 2009
Version: 2011-06-12
Modifications:
2009-09-13 (NM): First version of this code
2011-06-12 (NM): Added more flexible output file names
Notes: Input files must be evenly-sampled time series files in SAC
binary format (assuming NVHDR=6, IFTYPE=ITIME (1), LEVEN=1).
*************************************************************************/
#include "nmsaclib.h" /* SAC-related functions, and other stuff */
/* Self-documentation */
const char *sdoc[] = {
" ",
" SACBUTTER - SAC Butterworth high-pass and/or low-pass filter ",
" ",
" Usage: sacbutter [parameters] -f sac_files ",
" sacbutter <stdin [parameters] > stdout ",
" ",
" Parameters and defaults: ",
" -f list of SAC binary waveform files (or read from stdin) ",
" ",
" -b1 0 Butterworth low-cut frequency in Hz (0 = no filter) ",
" -b2 0 Butterworth high-cut frequency in Hz (0 = no filter) ",
" -bp 3 number of poles of Butterworth filters ",
" ",
" -bz flag: apply zerophase (two-pass) filter ",
" -o flag: overwrite input files (default appends \".flt\") ",
" ",
" -ext flt output file extension (if -f is set) ",
" -dir output directory (default is same as input directory) ",
" ",
" The parameters -flo and -fhi are acceptable aliases for -b1 and -b2, ",
" respectively. If the overwrite flag -o is set, original data are lost.",
" ",
" NM, 2011-06-12 ",
" ",
NULL};
/* Definitions */
#define RAD2DEG 57.29577951
int
main(int argc, char **argv)
{
/* Variables */
register int iarg, it; /* loop indices (files, time samples) */
SACHEAD hd; /* SAC header */
float *data=NULL; /* SAC input trace data */
char *outfname=NULL; /* output file name */
int isfile=0; /* file names flag */
int isvalid=0; /* internal SAC trace validation flag */
int swap; /* byte-swapping flag */
int overwrite=0; /* overwrite flag */
int zerophase=0; /* zero-phase flag */
float fmin, fmax; /* Butterworth filter frequencies */
int npoles; /* number of poles of Butterworth filters */
char *outdir=NULL; /* optional output directory name */
char *outext=NULL; /* output file extension */
/* Print documentation */
if (argc==1) {
int i=0;
while (sdoc[i]) fprintf(stderr, "%s\n", sdoc[i++]);
return (EXIT_FAILURE);
}
/* Command line parameters */
if (!fgetpar(argc, argv, "-flo", &fmin) && \
!fgetpar(argc, argv, "-b1", &fmin)) fmin=0.0;
if (!fgetpar(argc, argv, "-fhi", &fmax) && \
!fgetpar(argc, argv, "-b2", &fmax)) fmax=0.0;
if (!igetpar(argc, argv, "-bp", &npoles)) npoles=3;
if (!(zerophase=getflag(argc, argv, "-z")) && \
!(zerophase=getflag(argc, argv, "-bz"))) zerophase=0;
if (!(overwrite=getflag(argc, argv, "-o"))) overwrite=0;
if (!(isfile=getflag(argc, argv, "-f"))) isfile=0;
if (!sgetpar(argc, argv, "-ext", &outext)) outext="flt";
if (!sgetpar(argc, argv, "-dir", &outdir)) outdir=NULL;
/* Check parameters */
if (npoles<=0) error("%s: number of poles must be positive\n", argv[0]);
if (!isfile && isatty(1)) \
error("%s: can't write a SAC binary file to tty\n", argv[0]);
/* Loop over SAC files */
for (iarg=isfile+1; (iarg<argc && argv[iarg][0]!='-') || !isfile; iarg++) {
/* Read SAC binary file */
if (isfile) {
if (!read_sacbin_file(argv[iarg], &hd, &data, &swap)) \
error("%s: can't read %s\n", argv[0], argv[iarg]);
}
else {
if (!read_sacbin(stdin, &hd, &data, &swap)) break;
}
/* Check for valid trace type and sampling rate */
isvalid = (!hd.leven || hd.iftype==IRLIM || hd.iftype==IAMPH) ? 0 : 1;
if (hd.delta<=0.0) {
fprintf(stderr,"%s: invalid delta=%g, skipping %s\n", \
argv[0], hd.delta, (isfile) ? argv[iarg] : "stdin");
isvalid = 0;
}
if (isvalid) {
/* Subtract mean value from trace */
if (hd.depmen!=SAC_HEADER_FLOAT_UNDEFINED) {
for (it=0; it<hd.npts; it++) data[it] -= hd.depmen;
}
else {
remove_mean(data, hd.npts);
}
/* Butterworth low-cut and/or high-cut filter(s) */
if (npoles>0) {
butterworth_bandpass(npoles, fmin, fmax, hd.delta, hd.npts, \
zerophase, data);
}
/* Re-compute depmen/depmin/depmax */
set_sac_depminmax(&hd, data);
/* Write SAC binary files... */
if (isfile) {
if (overwrite) {
outfname = argv[iarg];
}
else {
if (!out_file_name(argv[iarg], outext, outdir, &outfname)) \
error("%s: can't set output file name\n", argv[0]);
}
if (!write_sacbin_file(outfname, hd, data, swap))\
error("%s: can't write %s\n", argv[0], outfname);
if (!overwrite) free(outfname);
}
/* ... or write SAC binary to stdout */
else {
if (!write_sacbin(stdout, hd, data, swap)) \
error("%s: can't write SAC trace to stdout\n", argv[0]);
}
}
}
return EXIT_SUCCESS;
}
/* END OF FILE */