Loading [MathJax]/extensions/TeX/AMSsymbols.js
LAL 7.7.0.1-da3b9d3
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
ButterworthTimeSeries.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Jolien Creighton, Teviet Creighton
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20#include <complex.h>
21#include <lal/LALStdlib.h>
22#include <lal/LALConstants.h>
23#include <lal/AVFactories.h>
24#include <math.h>
25#include <lal/IIRFilter.h>
26#include <lal/BandPassTimeSeries.h>
27
28/**
29 * \defgroup ButterworthTimeSeries_c Module ButterworthTimeSeries.c
30 * \ingroup BandPassTimeSeries_h
31 *
32 * \author Creighton, T. D.
33 *
34 * \brief Applies a low- or high-pass Butterworth filter to a time series.
35 *
36 * ### Synopsis ###
37 *
38 * The \ref ButterworthTimeSeries_c provides specific advice and
39 * guidelines for building a numerically stable time-domain filter,
40 * but the general procedure is as follows.
41 * <ol>
42 * <li>Decide on the desired filter response, and express it as a
43 * rational function of the frequency variable \f$w=\tan(\pi f\Delta
44 * t)\f$ (which maps the Nyquist interval onto the positive real
45 * axis).</li>
46 * <li>Factorize this rational function into zeros and poles,
47 * restricting oneself to the upper half of the \f$w\f$ complex plane.
48 * Assign these to one or more objects of type
49 * <tt><datatype>ZPGFilter</tt></li>
50 * <li>Use <tt>WToZ<datatype>ZPGFilter()</tt> to transform the
51 * filter to the more conventional \f$z=\exp(2\pi if\Delta t)\f$
52 * frequency variable.</li>
53 * <li>Use the routines in \ref IIRFilter_h to create IIR filters from
54 * the ZPG filters, and to apply them to data.</li>
55 * </ol>
56 *
57 * ### Description ###
58 *
59 * These routines perform an in-place time-domain band-pass filtering of
60 * a data sequence <tt>*series</tt>, using a Butterworth filter generated
61 * from parameters <tt>*params</tt>. The routines construct a filter with
62 * the square root of the desired amplitude response, which it then
63 * applied to the data once forward and once in reverse. This gives the
64 * full amplitude response with little or no frequency-dependent phase
65 * shift.
66 *
67 * The routine <tt>LALDButterworthREAL4TimeSeries()</tt> applies a
68 * double-precision filter to single-precision data, using
69 * <tt>LALDIIRFilterREAL4Vector()</tt> and
70 * <tt>LALDIIRFilterREAL4VectorR()</tt>.
71 *
72 * ### Algorithm ###
73 *
74 * The frequency response of a Butterworth low-pass filter is easiest to
75 * express in terms of the transformed frequency variable \f$w=\tan(\pi
76 * f\Delta t)\f$, where \f$\Delta t\f$ is the sampling interval (i.e.
77 * <tt>series->deltaT</tt>). In this parameter, then, the \e power
78 * response (attenuation) of the filter is:
79 * \f[
80 * |R|^2 = \sqrt{a} = \frac{1}{1+(w/w_c)^{2n}} \; ,
81 * \f]
82 * where \f$n\f$ is the filter order and \f$w_c\f$ is the characteristic
83 * frequency. We have written the attenuation as \f$\sqrt{a}\f$ to emphasize
84 * that the full attenuation \f$a\f$ is achieved only after filtering twice
85 * (once forward, once in reverse). Similarly, a Butterworth high-pass
86 * filter is given by
87 * \f[
88 * |R|^2 = \sqrt{a} = \frac{1}{1+(w_c/w)^{2n}} \; .
89 * \f]
90 * If one is given a filter order \f$n\f$, then the characteristic frequency
91 * can be determined from the attenuation at some any given frequency.
92 * Alternatively, \f$n\f$ and \f$w_c\f$ can both be computed given attenuations
93 * at two different frequencies.
94 *
95 * Frequencies in <tt>*params</tt> are assumed to be real frequencies \f$f\f$
96 * given in the inverse of the units used for the sampling interval
97 * <tt>series->deltaT</tt>. In order to be used, the pass band parameters
98 * must lie in the ranges given below; if a parameter lies outside of its
99 * range, then it is ignored and the filter is calculated from the
100 * remaining parameters. If too many parameters are missing, the routine
101 * will fail. The acceptable parameter ranges are:
102 *
103 * <dl>
104 * <dt>params->nMax</dt><dd> = 1, 2, \f$\ldots\f$</dd>
105 * <dt>params->f1, f2</dt><dd> \f$\in(0,\{2\times<tt>series->deltaT</tt>\}^{-1}) \f$</dd>
106 * <dt>params->a1, a2</dt><dd> \f$\in (0,1) \f$</dd>
107 * </dl>
108 *
109 * If both pairs of frequencies and amplitudes are given, then \c a1,
110 * \c a2 specify the minimal requirements on the attenuation of the
111 * filter at frequencies \c f1, \c f2. Whether the filter is a
112 * low- or high-pass filter is determined from the relative sizes of
113 * these parameters. In this case the \c nMax parameter is optional;
114 * if given, it specifies an upper limit on the filter order. If the
115 * desired attenuations would require a higher order, then the routine
116 * will sacrifice performance in the stop band in order to remain within
117 * the specified \c nMax.
118 *
119 * If one of the frequency/attenuation pairs is missing, then the filter
120 * is computed using the remaining pair and \c nMax (which must be
121 * given). The filter is taken to be a low-pass filter if \c f1,
122 * \c a1 are given, and high-pass if \c f2, \c a2 are given.
123 * If only one frequency and no corresponding attenuation is specified,
124 * then it is taken to be the characteristic frequency (i.e. the
125 * corresponding attenuation is assumed to be \f$\sqrt{a}=1/2\f$). If none
126 * of these conditions are met, the routine will return an error.
127 *
128 * The following table summarizes the decision algorithm. A \f$\bullet\f$
129 * symbol indicates that the parameter is specified in the range given
130 * above. A \f$\circ\f$ symbol indicates that the parameter is <i>not given</i>, i.e. not specified in the valid range.
131 *
132 * <table>
133 * <tr><th>nMax</th><th>f1</th><th>a1</th><th>f2</th><th>a2</th><th>Procedure</th></tr>
134 * <tr><td>\f$\circ\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>
135 * Type of filter (low- or high-pass), \f$w_c\f$, and \f$n\f$ are
136 * computed from all four transition-band parameters.</td></tr>
137 * <tr><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>
138 * Ditto, but if the resulting \f$n>\f$ \c nMax, \f$w_c\f$ is
139 * computed from \c nMax and the (\c f,\c a)
140 * pair with the \e larger \c a.</td></tr>
141 * <tr><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>\f$\circ\f$</td><td>
142 * Low-pass filter; \f$w_c\f$ is computed from \c nMax,
143 * \c f1, and \c a1.</td></tr>
144 * <tr><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>\f$\bullet\f$</td><td>
145 * Ditto; \c a2 is ignored.</td></tr>
146 * <tr><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>
147 * Ditto; \c f2 is ignored.</td></tr>
148 * <tr><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>\f$\circ\f$</td><td>\f$\circ\f$</td><td>
149 * Low-pass filter; \f$w_c\f$ is computed as above with \c a1
150 * treated as 1/4.</td></tr>
151 * <tr><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>\f$\circ\f$</td><td>\f$\bullet\f$</td><td>
152 * Ditto; \c a2 is ignored.</td></tr>
153 * <tr><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>\f$\circ\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>
154 * High-pass filter; \f$w_c\f$ is computed from \c nMax,
155 * \c f2, and \c a2.</td></tr>
156 * <tr><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>
157 * Ditto; \c a1 is ignored.</td></tr>
158 * <tr><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>
159 * Ditto; \c f1 is ignored.</td></tr>
160 * <tr><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>\f$\circ\f$</td><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>
161 * High-pass filter; \f$w_c\f$ is computed as above with \c a2
162 * treated as 1/4.</td></tr>
163 * <tr><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>\f$\bullet\f$</td><td>\f$\bullet\f$</td><td>\f$\circ\f$</td><td>
164 * Ditto; \c a1 is ignored.</td></tr>
165 * <tr><td colspan=5>Other</td><td>Subroutine returns an error.</td></tr>
166 * </table>
167 *
168 * Once an order \f$n\f$ and characteristic frequency \f$w_c\f$ are known, the
169 * zeros and poles of a ZPG filter are readily determined. A stable,
170 * physically realizable Butterworth filter will have \f$n\f$ poles evenly
171 * spaced on the upper half of a circle of radius \f$w_c\f$; that is,
172 * \f[
173 * R = \frac{(-iw_c)^n}{\prod_{k=0}^{n-1}(w - w_c e^{2\pi i(k+1/2)/n})}
174 * \f]
175 * for a low-pass filter, and
176 * \f[
177 * R = \frac{w^n}{\prod_{k=0}^{n-1}(w - w_c e^{2\pi i(k+1/2)/n})}
178 * \f]
179 * for a high-pass filter. By choosing only poles on the upper-half
180 * plane, one ensures that after transforming to \f$z\f$ the poles will have
181 * \f$|z|<1\f$. Furthermore, the phase factor \f$(-i)^n\f$ in the numerator of
182 * the low-pass filter is chosen so that the DC response is purely real;
183 * this ensures that the response function in the \f$z\f$-plane will have a
184 * real gain factor, and the resulting IIR filter will be physically
185 * realizable. The high-pass filter has a purely real response at
186 * Nyquist (\f$w\rightarrow\infty\f$), which similarly gives a physical IIR
187 * filter.
188 *
189 * Although higher orders \f$n\f$ would appear to produce better (i.e.
190 * sharper) filter responses, one rapidly runs into numerical errors, as
191 * one ends up adding and subtracting \f$n\f$ large numbers to obtain small
192 * filter responses. One way around this is to break the filter up into
193 * several lower-order filters. The routines in this module do just
194 * that. Poles are paired up across the imaginary axis, (and combined
195 * with pairs of zeros at \f$w=0\f$ for high-pass filters,) to form \f$[n/2]\f$
196 * second-order filters. If \f$n\f$ is odd, there will be an additional
197 * first-order filter, with one pole at \f$w=iw_c\f$ (and one zero at \f$w=0\f$
198 * for a high-pass filter).
199 *
200 * Each ZPG filter in the \f$w\f$-plane is first transformed to the \f$z\f$-plane
201 * by a bilinear transformation, and is then used to construct a
202 * time-domain IIR filter. Each filter is then applied to the time
203 * series. As mentioned in the description above, the filters are
204 * designed to give an overall amplitude response that is the square root
205 * of the desired attenuation; however, each time-domain filter is
206 * applied to the data stream twice: once in the normal sense, and once
207 * in the time-reversed sense. This gives the full attenuation with very
208 * little frequency-dependent phase shift.
209 *
210 */
211/** @{ */
212
213/* Prototype for a local input parsing routine. */
214static int
216 INT4 *n,
217 REAL8 *wc,
218 REAL8 deltaT );
219
220
221#undef COMPLEX_DATA
222#undef SINGLE_PRECISION
223
224#define COMPLEX_DATA
225#define SINGLE_PRECISION
227#undef SINGLE_PRECISION
229#undef COMPLEX_DATA
230#define SINGLE_PRECISION
232#undef SINGLE_PRECISION
234
235/**
236 * Deprecated.
237 * \deprecated Use XLALButterworthREAL4TimeSeries() instead.
238 */
239void
241 REAL4TimeSeries *series,
242 PassBandParamStruc *params )
243{
244 INT4 n; /* The filter order. */
245 INT4 type; /* The pass-band type: high, low, or undeterminable. */
246 INT4 i; /* An index. */
247 INT4 j; /* Another index. */
248 REAL8 wc; /* The filter's transformed frequency. */
249
250 INITSTATUS(stat);
251 ATTATCHSTATUSPTR(stat);
252
253 /* Make sure the input pointers are non-null. */
255 BANDPASSTIMESERIESH_MSGENUL);
257 BANDPASSTIMESERIESH_MSGENUL);
259 BANDPASSTIMESERIESH_MSGENUL);
261 BANDPASSTIMESERIESH_MSGENUL);
262
263 /* Parse the pass-band parameter structure. I separate this into a
264 local static subroutine because it's an icky mess of conditionals
265 that would clutter the logic of the main routine. */
266 type=XLALParsePassBandParamStruc(params,&n,&wc,series->deltaT);
267 if(type<0) {
269 ABORT(stat,BANDPASSTIMESERIESH_EBAD,BANDPASSTIMESERIESH_MSGEBAD);
270 }
271
272 /* An order n Butterworth filter has n poles spaced evenly along a
273 semicircle in the upper complex w-plane. By pairing up poles
274 symmetric across the imaginary axis, the filter gan be decomposed
275 into [n/2] filters of order 2, plus perhaps an additional order 1
276 filter. The following loop pairs up poles and applies the
277 filters with order 2. */
278 for(i=0,j=n-1;i<j;i++,j--){
279 REAL4 theta=LAL_PI*(i+0.5)/n;
280 REAL4 ar=wc*cos(theta);
281 REAL4 ai=wc*sin(theta);
282 REAL4IIRFilter *iirFilter=NULL;
283 COMPLEX8ZPGFilter *zpgFilter=NULL;
284
285 /* Generate the filter in the w-plane. */
286 if(type==2){
287 TRY(LALCreateCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter,2,2),
288 stat);
289 zpgFilter->zeros->data[0]=0.0;
290 zpgFilter->zeros->data[1]=0.0;
291 zpgFilter->gain=1.0;
292 }else{
293 TRY(LALCreateCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter,0,2),
294 stat);
295 zpgFilter->gain=-wc*wc;
296 }
297 zpgFilter->poles->data[0]=ar;
298 zpgFilter->poles->data[0]+=ai*I;
299 zpgFilter->poles->data[1]=-ar;
300 zpgFilter->poles->data[1]+=ai*I;
301
302 /* Transform to the z-plane and create the IIR filter. */
303 LALWToZCOMPLEX8ZPGFilter(stat->statusPtr,zpgFilter);
304 BEGINFAIL(stat)
305 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),
306 stat);
307 ENDFAIL(stat);
308 LALCreateREAL4IIRFilter(stat->statusPtr,&iirFilter,zpgFilter);
309 BEGINFAIL(stat)
310 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),
311 stat);
312 ENDFAIL(stat);
313
314 /* Filter the data, once each way. */
315 LALIIRFilterREAL4Vector(stat->statusPtr,series->data,iirFilter);
316 BEGINFAIL(stat) {
317 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),
318 stat);
319 TRY(LALDestroyREAL4IIRFilter(stat->statusPtr,&iirFilter),stat);
320 } ENDFAIL(stat);
321 LALIIRFilterREAL4VectorR(stat->statusPtr,series->data,iirFilter);
322 BEGINFAIL(stat) {
323 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),
324 stat);
325 TRY(LALDestroyREAL4IIRFilter(stat->statusPtr,&iirFilter),stat);
326 } ENDFAIL(stat);
327
328 /* Free the filters. */
329 TRY(LALDestroyREAL4IIRFilter(stat->statusPtr,&iirFilter),stat);
330 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),stat);
331 }
332
333 /* Next, this conditional applies the possible order 1 filter
334 corresponding to an unpaired pole on the imaginary w axis. */
335 if(i==j){
336 REAL4IIRFilter *iirFilter=NULL;
337 COMPLEX8ZPGFilter *zpgFilter=NULL;
338
339 /* Generate the filter in the w-plane. */
340 if(type==2){
341 TRY(LALCreateCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter,1,1),
342 stat);
343 *zpgFilter->zeros->data=0.0;
344 zpgFilter->gain=1.0;
345 }else{
346 TRY(LALCreateCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter,0,1),
347 stat);
348 zpgFilter->gain=-wc*I;
349 }
350 *zpgFilter->poles->data=wc*I;
351
352 /* Transform to the z-plane and create the IIR filter. */
353 LALWToZCOMPLEX8ZPGFilter(stat->statusPtr,zpgFilter);
354 BEGINFAIL(stat)
355 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),
356 stat);
357 ENDFAIL(stat);
358 LALCreateREAL4IIRFilter(stat->statusPtr,&iirFilter,zpgFilter);
359 BEGINFAIL(stat)
360 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),
361 stat);
362 ENDFAIL(stat);
363
364 /* Filter the data, once each way. */
365 LALIIRFilterREAL4Vector(stat->statusPtr,series->data,iirFilter);
366 BEGINFAIL(stat) {
367 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),
368 stat);
369 TRY(LALDestroyREAL4IIRFilter(stat->statusPtr,&iirFilter),stat);
370 } ENDFAIL(stat);
371 LALIIRFilterREAL4VectorR(stat->statusPtr,series->data,iirFilter);
372 BEGINFAIL(stat) {
373 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),
374 stat);
375 TRY(LALDestroyREAL4IIRFilter(stat->statusPtr,&iirFilter),stat);
376 } ENDFAIL(stat);
377
378 /* Free the filters. */
379 TRY(LALDestroyREAL4IIRFilter(stat->statusPtr,&iirFilter),stat);
380 TRY(LALDestroyCOMPLEX8ZPGFilter(stat->statusPtr,&zpgFilter),stat);
381 }
382
383 /* Normal exit. */
384 DETATCHSTATUSPTR(stat);
385 RETURN(stat);
386}
387
388
389/**
390 * Deprecated.
391 * \deprecated Use XLALButterworthREAL8TimeSeries() instead.
392 */
393void
395 REAL8TimeSeries *series,
396 PassBandParamStruc *params )
397{
398 INITSTATUS(stat);
399
400 if (XLALButterworthREAL8TimeSeries(series,params)<0)
401 {
402 int code=xlalErrno;
404 switch(code){
405 case XLAL_EFAULT:
406 ABORT(stat,BANDPASSTIMESERIESH_ENUL,BANDPASSTIMESERIESH_MSGENUL);
407 case XLAL_EINVAL:
408 ABORT(stat,BANDPASSTIMESERIESH_EBAD,BANDPASSTIMESERIESH_MSGEBAD);
409 default:
410 ABORTXLAL(stat);
411 }
412 }
413
414 RETURN(stat);
415}
416
417/**
418 * Deprecated.
419 * \deprecated
420 */
421void
423 REAL4TimeSeries *series,
424 PassBandParamStruc *params )
425{
426 INITSTATUS(stat);
427
428 if (XLALButterworthREAL4TimeSeries(series,params)<0)
429 {
430 int code=xlalErrno;
432 switch(code){
433 case XLAL_EFAULT:
434 ABORT(stat,BANDPASSTIMESERIESH_ENUL,BANDPASSTIMESERIESH_MSGENUL);
435 case XLAL_EINVAL:
436 ABORT(stat,BANDPASSTIMESERIESH_EBAD,BANDPASSTIMESERIESH_MSGEBAD);
437 default:
438 ABORTXLAL(stat);
439 }
440 }
441
442 RETURN(stat);
443}
444
445
446static int
448 INT4 *n,
449 REAL8 *wc,
450 REAL8 deltaT )
451 /* This local function parses the pass band parameters according
452 to the rules given in the documentation to this module,
453 computing the order and characteristic frequency (in the
454 w-variable) of the desired filter. The function returns 2 for
455 a high-pass filter, 1 for a low-pass filter, and -1 if the
456 parameters were poorly specified. */
457{
458 /* Okay, first define all the temporary variables, and figure out
459 just which parameter combinations have been given. */
460 REAL8 w1=deltaT*params->f1; /* Dimensionless frequencies; */
461 REAL8 w2=deltaT*params->f2; /* later transformed to w-plane. */
462 REAL8 a1=params->a1; /* Attenuation factors; later a */
463 REAL8 a2=params->a2; /* square root will be taken. */
464 BOOLEAN w1Given = (w1>0.0)&&(w1<0.5);
465 BOOLEAN w2Given = (w2>0.0)&&(w2<0.5);
466 BOOLEAN both1Given = w1Given && (a1>0.0)&&(a1<1.0);
467 BOOLEAN both2Given = w2Given && (a2>0.0)&&(a2<1.0);
468
469 /* If both frequencies and both attenuations have been given,
470 compute the required filter order and characteristic frequency
471 from them. */
472 if(both1Given && both2Given){
473 REAL8 aLow = w1<w2 ? sqrt(a1) : sqrt(a2);
474 REAL8 aHigh = w1>w2 ? sqrt(a1) : sqrt(a2);
475 REAL8 wLow = w1<w2 ? tan(LAL_PI*w1) : tan(LAL_PI*w2);
476 REAL8 wHigh = w1>w2 ? tan(LAL_PI*w1) : tan(LAL_PI*w2);
477
478 /* First make sure that two different frequencies and attenuations
479 have been specified. */
480 if(w1==w2){
481 XLALPrintInfo("XLAL Info - %s: ", __func__ );
482 XLALPrintInfo("The two frequencies defining the transition band"
483 " are the same.\n");
485 }
486 if(a1==a2){
487 XLALPrintInfo("XLAL Info - %s: ", __func__ );
488 XLALPrintInfo("The two attenuations across the transition band"
489 " are the same.\n");
491 }
492
493 /* Now compute the filter order. */
494 if(aHigh>aLow) /* High-pass. */
495 *n=(INT4)(0.5*log((1.0/aHigh-1.0)/(1.0/aLow-1.0))
496 /log(wLow/wHigh))+1;
497 else /* Low-pass. */
498 *n=(INT4)(0.5*log((1.0/aHigh-1.0)/(1.0/aLow-1.0))
499 /log(wHigh/wLow))+1;
500
501 /* If a positive params->nMax less than *n has been specified,
502 reduce to that order, with appropriate warnings. */
503 if((params->nMax>0)&&(params->nMax<*n)){
504 XLALPrintWarning("XLAL Warning - %s: ", __func__ );
505 XLALPrintWarning("Filter order required to achieve requested"
506 " performance exceeds\n"
507 "\tspecified limit.\n");
508 XLALPrintWarning("\tRequired: %i Limit: %i\n",*n,params->nMax);
509 *n=params->nMax;
510 }
511
512 /* Compute the characteristic frequency of the filter using the
513 filter order and the minimal attenuation in the pass-band. */
514 if(aHigh>aLow){ /* High-pass */
515 *wc=wHigh*pow(1.0/aHigh - 1.0,-0.5/((REAL8)(*n)));
516 return 2;
517 }else{ /* Low-pass */
518 *wc=wLow*pow(1.0/aLow - 1.0,0.5/((REAL8)(*n)));
519 return 1;
520 }
521 }
522
523 /* If one end of the transition band has both frequency and
524 attenuation specified, and the other does not, use only the end
525 that does. The filter order must also be specified for this and
526 all future cases. */
527 else{
528 if(params->nMax<=0){
529 XLALPrintInfo("XLAL Info - %s: ", __func__ );
530 XLALPrintInfo("The filter order must be given if one or more"
531 " frequencies or\n"
532 "\tattenuations of the transition band have not been"
533 " (validly) specified.");
535 }
536 if(both1Given){
537 REAL8 wLow=tan(LAL_PI*w1);
538 REAL8 aLow=sqrt(a1);
539 *n=params->nMax;
540 *wc=wLow*pow(1.0/aLow - 1.0, -0.5/((REAL8)(*n)));
541 return 1;
542 }else if(both2Given){
543 REAL8 wHigh=tan(LAL_PI*w2);
544 REAL8 aHigh=sqrt(a2);
545 *n=params->nMax;
546 *wc=wHigh*pow(1.0/aHigh - 1.0, 0.5/((REAL8)(*n)));
547 return 2;
548 }
549
550 /* If no attenuations are given, then there had better be only one
551 frequency given, otherwise we don't know whether to make a low-
552 or a high-pass filter. */
553 else if(w1Given && w2Given){
554 XLALPrintInfo("XLAL Info - %s: ", __func__ );
555 XLALPrintInfo("Neither attenuation has been specified, so only"
556 " one frequency should be.");
558 }
559
560 /* Treat the given frequency as the characteristic frequency and
561 the specified nMax as the filter order. */
562 else{
563 *n=params->nMax;
564 if(w2Given){
565 *wc=tan(LAL_PI*w2);
566 return 2;
567 }else if(w1Given){
568 *wc=tan(LAL_PI*w1);
569 return 1;
570 }else{
571 XLALPrintInfo("XLAL Info - %s: ", __func__ );
572 XLALPrintInfo("No frequencies within the Nyquist band have"
573 " been specified!");
575 }
576 }
577 }
578}
579/** @} */
int XLALButterworthREAL8TimeSeries(REAL8TimeSeries *series, PassBandParamStruc *params)
int XLALButterworthREAL4TimeSeries(REAL4TimeSeries *series, PassBandParamStruc *params)
#define ABORT(statusptr, code, mesg)
#define ENDFAIL(statusptr)
#define TRY(func, statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define ASSERT(assertion, statusptr, code, mesg)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define ABORTXLAL(sp)
#define BEGINFAIL(statusptr)
int XLALPrintWarning(const char *fmt,...)
Definition: XLALError.c:79
int XLALPrintInfo(const char *fmt,...)
Definition: XLALError.c:90
#define BANDPASSTIMESERIESH_ENUL
Unexpected null pointer in arguments.
#define BANDPASSTIMESERIESH_EBAD
Bad filter parameters.
void LALWToZCOMPLEX8ZPGFilter(LALStatus *stat, COMPLEX8ZPGFilter *filter)
Deprecated.
void LALButterworthREAL8TimeSeries(LALStatus *stat, REAL8TimeSeries *series, PassBandParamStruc *params)
Deprecated.
void LALDButterworthREAL4TimeSeries(LALStatus *stat, REAL4TimeSeries *series, PassBandParamStruc *params)
Deprecated.
void LALButterworthREAL4TimeSeries(LALStatus *stat, REAL4TimeSeries *series, PassBandParamStruc *params)
Deprecated.
static int XLALParsePassBandParamStruc(PassBandParamStruc *params, INT4 *n, REAL8 *wc, REAL8 deltaT)
void LALCreateREAL4IIRFilter(LALStatus *stat, REAL4IIRFilter **output, COMPLEX8ZPGFilter *input)
Deprecated.
void LALCreateCOMPLEX8ZPGFilter(LALStatus *stat, COMPLEX8ZPGFilter **output, INT4 numZeros, INT4 numPoles)
void LALDestroyREAL4IIRFilter(LALStatus *stat, REAL4IIRFilter **input)
Deprecated.
void LALDestroyCOMPLEX8ZPGFilter(LALStatus *stat, COMPLEX8ZPGFilter **input)
void LALIIRFilterREAL4Vector(LALStatus *status, REAL4Vector *vector, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
void LALIIRFilterREAL4VectorR(LALStatus *status, REAL4Vector *vector, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
#define LAL_PI
Archimedes's constant, pi.
Definition: LALConstants.h:179
unsigned char BOOLEAN
Boolean logical type, see Headers LAL(Atomic)Datatypes.h for more details.
double REAL8
Double precision real floating-point number (8 bytes).
int32_t INT4
Four-byte signed integer.
float REAL4
Single precision real floating-point number (4 bytes).
#define xlalErrno
Modifiable lvalue containing the XLAL error number.
Definition: XLALError.h:571
#define XLAL_ERROR(...)
Macro to invoke a failure from a XLAL routine returning an integer.
Definition: XLALError.h:700
int XLALClearErrno(void)
Clears the XLAL error number, returns the old value.
Definition: XLALError.c:363
@ XLAL_EFAULT
Invalid pointer.
Definition: XLALError.h:408
@ XLAL_EINVAL
Invalid argument.
Definition: XLALError.h:409
COMPLEX8 * data
Pointer to the data array.
Definition: LALDatatypes.h:168
See DATATYPE-ZPGFilter types for details.
Definition: LALDatatypes.h:921
COMPLEX8Vector * poles
Definition: LALDatatypes.h:925
COMPLEX8Vector * zeros
Definition: LALDatatypes.h:924
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
struct tagLALStatus * statusPtr
Pointer to the next node in the list; NULL if this function is not reporting a subroutine error.
Definition: LALDatatypes.h:954
This structure stores data used for constructing a low- or high-pass filter: either the order and cha...
REAL8 a1
The minimal desired attenuation factors at the reference frequencies.
REAL8 a2
The minimal desired attenuation factors at the reference frequencies.
REAL8 f1
The reference frequencies of the transition band.
INT4 nMax
The maximum desired filter order (actual order may be less if specified attenuations do not require a...
REAL8 f2
The reference frequencies of the transition band.
This structure stores the direct and recursive REAL4 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:152
Time series of REAL4 data, see DATATYPE-TimeSeries types for more details.
Definition: LALDatatypes.h:570
REAL4Sequence * data
The sequence of sampled data.
Definition: LALDatatypes.h:576
REAL8 deltaT
The time step between samples of the time series in seconds.
Definition: LALDatatypes.h:573
REAL4 * data
Pointer to the data array.
Definition: LALDatatypes.h:150
Time series of REAL8 data, see DATATYPE-TimeSeries types for more details.
Definition: LALDatatypes.h:580