LAL  7.5.0.1-8083555
IIRFilter.h
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 #ifndef _IIRFILTER_H
21 #define _IIRFILTER_H
22 
23 #include <lal/LALStdlib.h>
24 #include <lal/ZPGFilter.h>
25 
26 #if defined(__cplusplus)
27 extern "C" {
28 #elif 0
29 } /* so that editors will match preceding brace */
30 #endif
31 
32 /**
33  * \defgroup IIRFilter_h Header IIRFilter.h
34  * \ingroup lal_tdfilter
35  * \author Creighton, T. D.
36  *
37  * \brief Provides routines to make and apply IIR filters.
38  *
39  * ### Synopsis ###
40  *
41  * \code
42  * #include <lal/IIRFilter.h>
43  * \endcode
44  *
45  * The \ref IIRFilter_h provides routines for creating actual
46  * time-domain filters from the ZPG representation, and applying these
47  * filters to data.
48  *
49  * This header covers routines that create, destroy, and apply
50  * generic time-domain filters, given by objects of type
51  * <tt><datatype>IIRFilter</tt>, where <tt><datatype></tt> is either
52  * \c REAL4 or \c REAL8.
53  *
54  * An IIR (Infinite Impulse Response) filter is a generalized linear
55  * causal time-domain filter, in which the filter output \f$y_n=y(t_n)\f$ at
56  * any sampled time \f$t_n=t_0+n\Delta t\f$ is a linear combination of the
57  * input \f$x\f$ \e and output \f$y\f$ at previous sampled times:
58  * \f[
59  * y_n = \sum_{k=0}^M c_k x_{n-k} + \sum_{l=1}^N d_l y_{n-l} \; .
60  * \f]
61  * The coefficients \f$c_k\f$ are called the direct filter coefficients, and
62  * the coefficients \f$d_l\f$ are the recursive filter coefficients. The
63  * filter order is the larger of \f$M\f$ or \f$N\f$, and determines how far back
64  * in time the filter must look to determine its next output. However,
65  * the recursive nature of the filter means that the output can depend on
66  * input arbitrarily far in the past; hence the name "infinite impulse
67  * response". Nonetheless, for a well-designed, stable filter, the
68  * actual filter response to an impulse should diminish rapidly beyond
69  * some characteristic timescale.
70  *
71  * Note that nonrecursive FIR (Finite Impulse Response) filters are
72  * considered a subset of IIR filters, having \f$N=0\f$.
73  *
74  * For practical implementation, it is convenient to express the bilinear
75  * equation above as two linear equations involving an auxiliary sequence
76  * \f$w\f$:
77  * \f[
78  * w_n = x_n + \sum_{l=1}^N d_l w_{n-l} \; ,
79  * \f]
80  * \f[
81  * y_n = \sum_{k=0}^M c_k w_{n-k} \; .
82  * \f]
83  * The equivalence of this to the first expression is not obvious, but
84  * can be proven by mathematical induction. The advantage of the
85  * auxiliary variable representation is twofold. First, when one is
86  * feeding data point by point to the filter, the filter needs only
87  * "remember" the previous \f$M\f$ or \f$N\f$ (whichever is larger) values of
88  * \f$w\f$, rather than remembering the previous \f$M\f$ values of \f$x\f$ \e and
89  * the previous \f$N\f$ values of \f$y\f$. Second, when filtering a large stored
90  * data vector, the filter response can be computed in place: one first
91  * runs forward through the vector replacing \f$x\f$ with \f$w\f$, and then
92  * backward replacing \f$w\f$ with \f$y\f$.
93  *
94  * Although the IIR filters in these routines are explicitly real, one
95  * can consider formally their complex response. A sinusoidal input can
96  * thus be written as \f$x_n=X\exp(2\pi ifn\Delta t)=Xz^n\f$, where \f$X\f$ is a
97  * complex amplitude and \f$z=\exp(2\pi if\Delta t)\f$ is a complex
98  * parametrization of the frequency. By linearity, the output must also
99  * be sinusoidal: \f$y_m=Y\exp(2\pi ifm\Delta t)=Yz^m\f$. Putting these into
100  * the bilinear equation, one can easily compute the filter's complex
101  * transfer function:
102  * \f[
103  * T(z) = \frac{Y}{X} = \frac{\sum_{k=0}^M c_k z^{-k}}
104  * {1 - \sum_{l=1}^N d_l z^{-l}}
105  * \f]
106  * This can be readily converted to and from the "zeros, poles, gain"
107  * representation of a filter, which expresses \f$T(z)\f$ as a factored
108  * rational function of \f$z\f$.
109  *
110  * It should also be noted that, in the routines covered by this header,
111  * I have adopted the convention of including a redundant recursive
112  * coefficient \f$d_0\f$, in order to make the indexing more intuitive. For
113  * formal correctness \f$d_0\f$ should be set to \f$-1\f$, although the filtering
114  * routines never actually use this coefficient.
115  *
116  */
117 /** @{ */
118 
119 /**
120  * @{
121  * \defgroup CreateIIRFilter_c Module CreateIIRFilter.c
122  * \defgroup DestroyIIRFilter_c Module DestroyIIRFilter.c
123  * \defgroup IIRFilter_c Module IIRFilter.c
124  * \defgroup IIRFilterVector_c Module IIRFilterVector.c
125  * \defgroup IIRFilterVectorR_c Module IIRFilterVectorR.c
126  * @}
127  */
128 
129 /** \name Error Codes */
130 /** @{ */
131 #define IIRFILTERH_ENUL 1 /**< Unexpected null pointer in arguments */
132 #define IIRFILTERH_EOUT 2 /**< Output handle points to a non-null pointer */
133 #define IIRFILTERH_EMEM 3 /**< Memory allocation error */
134 #define IIRFILTERH_EPAIR 4 /**< Input has unpaired nonreal poles or zeros */
135 /** @} */
136 
137 /** \cond DONT_DOXYGEN */
138 #define IIRFILTERH_MSGENUL "Unexpected null pointer in arguments"
139 #define IIRFILTERH_MSGEOUT "Output handle points to a non-null pointer"
140 #define IIRFILTERH_MSGEMEM "Memory allocation error"
141 #define IIRFILTERH_MSGEPAIR "Input has unpaired nonreal poles or zeros"
142 /** \endcond */
143 
144 /**
145  * This structure stores the direct and recursive REAL4 filter coefficients, as
146  * well as the history of the auxiliary sequence \f$w\f$.
147  * The length of the history vector gives the order of the filter
148  */
149 #ifdef SWIG /* SWIG interface directives */
150 SWIGLAL(IMMUTABLE_MEMBERS(tagREAL4IIRFilter, name));
151 #endif /* SWIG */
152 typedef struct tagREAL4IIRFilter{
153  const CHAR *name; /**< User assigned name. */
154  REAL8 deltaT; /**< Sampling time interval of the filter; If \f$\leq0\f$, it will be ignored (ie it will be taken from the data stream) */
155  REAL4Vector *directCoef; /**< The direct filter coefficients. */
156  REAL4Vector *recursCoef; /**< The recursive filter coefficients. */
157  REAL4Vector *history; /**< The previous values of w. */
159 
160 /**
161  * This structure stores the direct and recursive REAL8 filter coefficients, as
162  * well as the history of the auxiliary sequence \f$w\f$.
163  * The length of the history vector gives the order of the filter
164  */
165 #ifdef SWIG /* SWIG interface directives */
166 SWIGLAL(IMMUTABLE_MEMBERS(tagREAL8IIRFilter, name));
167 #endif /* SWIG */
168 typedef struct tagREAL8IIRFilter{
169  const CHAR *name; /**< User assigned name. */
170  REAL8 deltaT; /**< Sampling time interval of the filter; If \f$\leq0\f$, it will be ignored (ie it will be taken from the data stream). */
171  REAL8Vector *directCoef; /**< The direct filter coefficients. */
172  REAL8Vector *recursCoef; /**< The recursive filter coefficients. */
173  REAL8Vector *history; /**< The previous values of w. */
175 
176 /**
177  * This structure stores the direct and recursive REAL8 filter coefficients, as
178  * well as the complex-valued history of the auxiliary sequence \f$w\f$.
179  * The length of the history vector gives the order of the filter
180  */
181 #ifdef SWIG /* SWIG interface directives */
182 SWIGLAL(IMMUTABLE_MEMBERS(tagCOMPLEX16IIRFilter, name));
183 #endif /* SWIG */
184 typedef struct tagCOMPLEX16IIRFilter{
185  const CHAR *name; /**< User assigned name. */
186  REAL8 deltaT; /**< Sampling time interval of the filter; If \f$\leq0\f$, it will be ignored (ie it will be taken from the data stream). */
187  REAL8Vector *directCoef; /**< The direct filter coefficients. */
188  REAL8Vector *recursCoef; /**< The recursive filter coefficients. */
189  COMPLEX16Vector *history; /**< The previous values of w. */
191 
192 /** @} */
193 
194 /* Function prototypes. */
201 
210 
213 /* WARNING: THIS FUNCTION IS OBSOLETE */
215 /* REAL8 LALDIIRFilter( REAL8 x, REAL8IIRFilter *filter ); */
216 #define LALDIIRFilter(x,f) XLALIIRFilterREAL8(x,f)
217 
218 
219 
220 /* ----- CreateIIRFilter.c ---------- */
221 void
224  COMPLEX8ZPGFilter *input );
225 
226 void
229  COMPLEX16ZPGFilter *input );
230 
231 /* ----- DestroyIIRFilter.c ---------- */
232 void
234  REAL4IIRFilter **input );
235 
236 void
238  REAL8IIRFilter **input );
239 
240 /* ----- IIRFilter.c ---------- */
241 void
243  REAL4 *output,
244  REAL4 input,
245  REAL4IIRFilter *filter );
246 
247 void
249  REAL8 *output,
250  REAL8 input,
251  REAL8IIRFilter *filter );
252 
253 /* ----- IIRFilterVector.c ---------- */
254 void
256  REAL4Vector *vector,
257  REAL4IIRFilter *filter );
258 
259 void
261  REAL8Vector *vector,
262  REAL8IIRFilter *filter );
263 
264 void
266  REAL4Vector *vector,
267  REAL8IIRFilter *filter );
268 
269 /* ----- IIRFilterVectorR.c ---------- */
270 void
272  REAL4Vector *vector,
273  REAL4IIRFilter *filter );
274 
275 void
277  REAL8Vector *vector,
278  REAL8IIRFilter *filter );
279 
280 void
282  REAL4Vector *vector,
283  REAL8IIRFilter *filter );
284 
285 
286 #if 0
287 { /* so that editors will match succeeding brace */
288 #elif defined(__cplusplus)
289 }
290 #endif
291 
292 #endif /* _IIRFILTER_H */
int XLALIIRFilterCOMPLEX8Vector(COMPLEX8Vector *vector, COMPLEX16IIRFilter *filter)
int XLALIIRFilterCOMPLEX16Vector(COMPLEX16Vector *vector, COMPLEX16IIRFilter *filter)
int XLALIIRFilterReverseCOMPLEX8Vector(COMPLEX8Vector *vector, COMPLEX16IIRFilter *filter)
int XLALIIRFilterReverseCOMPLEX16Vector(COMPLEX16Vector *vector, COMPLEX16IIRFilter *filter)
int XLALIIRFilterReverseREAL8Vector(REAL8Vector *vector, REAL8IIRFilter *filter)
int XLALIIRFilterREAL8Vector(REAL8Vector *vector, REAL8IIRFilter *filter)
int XLALIIRFilterReverseREAL4Vector(REAL4Vector *vector, REAL8IIRFilter *filter)
int XLALIIRFilterREAL4Vector(REAL4Vector *vector, REAL8IIRFilter *filter)
const char *const name
type name
Definition: UserInput.c:193
REAL4IIRFilter * XLALCreateREAL4IIRFilter(COMPLEX8ZPGFilter *input)
void LALCreateREAL4IIRFilter(LALStatus *status, REAL4IIRFilter **output, COMPLEX8ZPGFilter *input)
Deprecated.
REAL8IIRFilter * XLALCreateREAL8IIRFilter(COMPLEX16ZPGFilter *input)
COMPLEX16IIRFilter * XLALCreateCOMPLEX16IIRFilter(COMPLEX16ZPGFilter *input)
void LALCreateREAL8IIRFilter(LALStatus *status, REAL8IIRFilter **output, COMPLEX16ZPGFilter *input)
Deprecated.
void LALDestroyREAL4IIRFilter(LALStatus *status, REAL4IIRFilter **input)
Deprecated.
void XLALDestroyREAL4IIRFilter(REAL4IIRFilter *filter)
void XLALDestroyCOMPLEX16IIRFilter(COMPLEX16IIRFilter *filter)
void XLALDestroyREAL8IIRFilter(REAL8IIRFilter *filter)
void LALDestroyREAL8IIRFilter(LALStatus *status, REAL8IIRFilter **input)
Deprecated.
void LALIIRFilterREAL8(LALStatus *status, REAL8 *output, REAL8 input, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
Definition: IIRFilter.c:112
REAL4 XLALIIRFilterREAL4(REAL4 x, REAL8IIRFilter *filter)
Definition: IIRFilter.c:220
REAL8 XLALIIRFilterREAL8(REAL8 x, REAL8IIRFilter *filter)
Definition: IIRFilter.c:182
REAL4 LALSIIRFilter(REAL4 x, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
Definition: IIRFilter.c:143
void LALIIRFilterREAL4(LALStatus *status, REAL4 *output, REAL4 input, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
Definition: IIRFilter.c:55
void LALDIIRFilterREAL4Vector(LALStatus *status, REAL4Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
void LALIIRFilterREAL8Vector(LALStatus *status, REAL8Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
void LALIIRFilterREAL4Vector(LALStatus *status, REAL4Vector *vector, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
void LALDIIRFilterREAL4VectorR(LALStatus *status, REAL4Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
void LALIIRFilterREAL4VectorR(LALStatus *status, REAL4Vector *vector, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
void LALIIRFilterREAL8VectorR(LALStatus *status, REAL8Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
double REAL8
Double precision real floating-point number (8 bytes).
char CHAR
One-byte signed integer, see Headers LAL(Atomic)Datatypes.h for more details.
float REAL4
Single precision real floating-point number (4 bytes).
This structure stores the direct and recursive REAL8 filter coefficients, as well as the complex-valu...
Definition: IIRFilter.h:184
const CHAR * name
User assigned name.
Definition: IIRFilter.h:185
COMPLEX16Vector * history
The previous values of w.
Definition: IIRFilter.h:189
REAL8Vector * directCoef
The direct filter coefficients.
Definition: IIRFilter.h:187
REAL8Vector * recursCoef
The recursive filter coefficients.
Definition: IIRFilter.h:188
REAL8 deltaT
Sampling time interval of the filter; If , it will be ignored (ie it will be taken from the data stre...
Definition: IIRFilter.h:186
Vector of type COMPLEX16, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:172
See DATATYPE-ZPGFilter types for details.
Definition: LALDatatypes.h:930
Vector of type COMPLEX8, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:163
See DATATYPE-ZPGFilter types for details.
Definition: LALDatatypes.h:921
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
This structure stores the direct and recursive REAL4 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:152
REAL4Vector * history
The previous values of w.
Definition: IIRFilter.h:157
REAL8 deltaT
Sampling time interval of the filter; If , it will be ignored (ie it will be taken from the data stre...
Definition: IIRFilter.h:154
const CHAR * name
User assigned name.
Definition: IIRFilter.h:153
REAL4Vector * recursCoef
The recursive filter coefficients.
Definition: IIRFilter.h:156
REAL4Vector * directCoef
The direct filter coefficients.
Definition: IIRFilter.h:155
Vector of type REAL4, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:145
This structure stores the direct and recursive REAL8 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:168
const CHAR * name
User assigned name.
Definition: IIRFilter.h:169
REAL8Vector * history
The previous values of w.
Definition: IIRFilter.h:173
REAL8Vector * recursCoef
The recursive filter coefficients.
Definition: IIRFilter.h:172
REAL8Vector * directCoef
The direct filter coefficients.
Definition: IIRFilter.h:171
REAL8 deltaT
Sampling time interval of the filter; If , it will be ignored (ie it will be taken from the data stre...
Definition: IIRFilter.h:170
Vector of type REAL8, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:154
void output(int gps_sec, int output_type)
Definition: tconvert.c:440