LAL  7.1.7.1-2d066e5
IIRFilterVector.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 <lal/LALStdlib.h>
21 #include <lal/IIRFilter.h>
22 
23 /**
24  * \addtogroup IIRFilterVector_c
25  * \author Creighton, T. D.
26  *
27  * \brief Applies an IIR filter to a data stream.
28  *
29  * ### Description ###
30  *
31  * These functions apply a generic time-domain filter given by an object
32  * <tt>*filter</tt> of type \c REAL4IIRFilter or \c REAL8IIRFilter
33  * to a list <tt>*vector</tt> of data representing a time series. This is
34  * done in place using the auxiliary data series formalism described in
35  * \ref IIRFilter.h, so as to accomodate potentially large data series.
36  * To filter a piece of a larger dataset, the calling routine may pass a
37  * vector structure whose data pointer and length fields specify a subset
38  * of a larger vector.
39  *
40  * The routine <tt>LALDIIRFilterREAL4Vector()</tt> applies a
41  * double-precision filter to single-precision data. It makes a single
42  * pass through the data, continuously updating the filter history at
43  * each step rather than storing the auxiliary array in-place. This
44  * reduces roundoff error by keeping \e all intermediate results to
45  * double-precision.
46  *
47  * ### Algorithm ###
48  *
49  * The implementation of <tt>LALDIIRFilterREAL4Vector()</tt> not only has
50  * lower truncation errors than <tt>LALIIRFilterREAL4Vector()</tt>, but
51  * also appears to be more computationally efficient, for reasons I have
52  * not yet determined; see the documentation for \ref IIRFilterTest.c.
53  * These combine to suggest that <tt>LALDIIRFilterREAL4Vector()</tt> is the
54  * better overall algorithm for filtering \c REAL4Vectors.
55  *
56  */
57 /** @{ */
58 
59 #undef COMPLEX_DATA
60 #undef SINGLE_PRECISION
61 
62 #define COMPLEX_DATA
63 #define SINGLE_PRECISION
64 #include "IIRFilterVector_source.c"
65 #undef SINGLE_PRECISION
66 #include "IIRFilterVector_source.c"
67 #undef COMPLEX_DATA
68 #define SINGLE_PRECISION
69 #include "IIRFilterVector_source.c"
70 #undef SINGLE_PRECISION
71 #include "IIRFilterVector_source.c"
72 
73 /**
74  * WARNING: THIS FUNCTION IS OBSOLETE.
75  * \deprecated
76  */
77 void
79  REAL4Vector *vector,
80  REAL4IIRFilter *filter )
81 {
82  INT4 i; /* Loop counter for data vector. */
83  INT4 j; /* Index for filter coeficients. */
84  INT4 k; /* Index for filter history. */
85  INT4 length; /* Length of vector. */
86  REAL4 *data; /* Vector data. */
87  REAL8 datum; /* Temporary working variable. */
88  INT4 directOrder; /* Number of direct filter coefficients. */
89  INT4 recursOrder; /* Number of recursive filter coefficients. */
90  INT4 numHist; /* The number of history data. */
91  REAL4 *directCoef; /* Direct filter coefficients. */
92  REAL4 *recursCoef; /* Recursive filter coefficients. */
93  REAL4 *history; /* Filter history. */
94  REAL4 *temp=NULL; /* Temporary storage for the filter history. */
95 
96  INITSTATUS(stat);
97 
98  /* Make sure all the structures have been initialized. */
99  ASSERT(vector,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
100  ASSERT(filter,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
101  ASSERT(vector->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
102  ASSERT(filter->directCoef,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
103  ASSERT(filter->recursCoef,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
104  ASSERT(filter->history,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
105  ASSERT(filter->directCoef->data,stat,IIRFILTERH_ENUL,
106  IIRFILTERH_MSGENUL);
107  ASSERT(filter->recursCoef->data,stat,IIRFILTERH_ENUL,
108  IIRFILTERH_MSGENUL);
109  ASSERT(filter->history->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
110  length=vector->length;
111  data=vector->data;
112  directOrder=filter->directCoef->length;
113  recursOrder=filter->recursCoef->length;
114  numHist=filter->history->length;
115  directCoef=filter->directCoef->data;
116  recursCoef=filter->recursCoef->data;
117  history=filter->history->data;
118  temp=(REAL4 *)LALMalloc(numHist*sizeof(REAL4));
119  if ( !temp ) {
120  ABORT(stat,IIRFILTERH_EMEM,IIRFILTERH_MSGEMEM);
121  }
122 
123  /* Compute the auxiliary data series. */
124  for(i=0;(i<recursOrder)&&(i<length);i++,data++){
125  datum=*data;
126  for(j=1;j<=i;j++)
127  datum+=data[-j]*recursCoef[j];
128  for(k=0;j<recursOrder;j++,k++)
129  datum+=history[k]*recursCoef[j];
130  *data=datum;
131  }
132  for(;i<length;i++,data++){
133  datum=*data;
134  for(j=1;j<recursOrder;j++)
135  datum+=data[-j]*recursCoef[j];
136  *data=datum;
137  }
138  data--;
139 
140  /* Store the last few auxiliary data to the temporary history. */
141  for(k=numHist-1;k>=length;k--)
142  temp[k]=history[k-length];
143  for(;k>=0;k--)
144  temp[k]=data[-k];
145 
146  /* Compute the output data series. */
147  for(;i>directOrder;i--,data--){
148  datum=*data*directCoef[0];
149  for(j=1;j<directOrder;j++)
150  datum+=data[-j]*directCoef[j];
151  *data=datum;
152  }
153  for(;i>0;i--,data--){
154  datum=*data*directCoef[0];
155  for(j=1;j<i;j++)
156  datum+=data[-j]*directCoef[j];
157  for(k=0;j<directOrder;j++,k++)
158  datum+=history[k]*directCoef[j];
159  *data=datum;
160  }
161 
162  /* Update the filter history from the temporary history. */
163  for(k=0;k<numHist;k++)
164  history[k]=temp[k];
165  LALFree(temp);
166 
167  /* Normal exit */
168  RETURN(stat);
169 }
170 
171 
172 /**
173  * WARNING: THIS FUNCTION IS OBSOLETE.
174  * \deprecated
175  */
176 void
178  REAL8Vector *vector,
179  REAL8IIRFilter *filter )
180 {
181  INITSTATUS(stat);
182 
183  /* Make sure all the structures have been initialized. */
184  ASSERT(vector,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
185  ASSERT(filter,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
186  ASSERT(vector->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
187  ASSERT(filter->directCoef,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
188  ASSERT(filter->recursCoef,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
189  ASSERT(filter->history,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
190  ASSERT(filter->directCoef->data,stat,IIRFILTERH_ENUL,
191  IIRFILTERH_MSGENUL);
192  ASSERT(filter->recursCoef->data,stat,IIRFILTERH_ENUL,
193  IIRFILTERH_MSGENUL);
194  ASSERT(filter->history->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
195 
196  if ( XLALIIRFilterREAL8Vector( vector, filter ) < 0 )
197  {
198  XLALClearErrno();
199  ABORT(stat,IIRFILTERH_EMEM,IIRFILTERH_MSGEMEM);
200  }
201 
202  /* Normal exit */
203  RETURN(stat);
204 }
205 
206 
207 /**
208  * WARNING: THIS FUNCTION IS OBSOLETE.
209  * \deprecated
210  */
211 void
213  REAL4Vector *vector,
214  REAL8IIRFilter *filter )
215 {
216  INITSTATUS(stat);
217 
218  /* Make sure all the structures have been initialized. */
219  ASSERT(vector,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
220  ASSERT(filter,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
221  ASSERT(vector->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
222  ASSERT(filter->directCoef,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
223  ASSERT(filter->recursCoef,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
224  ASSERT(filter->history,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
225  ASSERT(filter->directCoef->data,stat,IIRFILTERH_ENUL,
226  IIRFILTERH_MSGENUL);
227  ASSERT(filter->recursCoef->data,stat,IIRFILTERH_ENUL,
228  IIRFILTERH_MSGENUL);
229  ASSERT(filter->history->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
230 
231  if ( XLALIIRFilterREAL4Vector( vector, filter ) < 0 )
232  {
233  XLALClearErrno();
234  ABORT(stat,IIRFILTERH_EMEM,IIRFILTERH_MSGEMEM);
235  }
236 
237  /* Normal exit */
238  RETURN(stat);
239 }
240 /** @} */
int XLALIIRFilterREAL4Vector(REAL4Vector *vector, REAL8IIRFilter *filter)
This structure stores the direct and recursive REAL8 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:168
REAL8Vector * directCoef
The direct filter coefficients.
Definition: IIRFilter.h:171
REAL8Vector * recursCoef
The recursive filter coefficients.
Definition: IIRFilter.h:172
REAL8 * data
Pointer to the data array.
Definition: LALDatatypes.h:159
This structure stores the direct and recursive REAL4 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:152
REAL4Vector * directCoef
The direct filter coefficients.
Definition: IIRFilter.h:155
#define LALFree(p)
Definition: LALMalloc.h:96
REAL8Vector * history
The previous values of w.
Definition: IIRFilter.h:173
REAL4Vector * history
The previous values of w.
Definition: IIRFilter.h:157
float REAL4
Single precision real floating-point number (4 bytes).
#define IIRFILTERH_ENUL
Unexpected null pointer in arguments.
Definition: IIRFilter.h:131
Vector of type REAL4, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:145
#define RETURN(statusptr)
#define INITSTATUS(statusptr)
int XLALClearErrno(void)
Clears the XLAL error number, returns the old value.
Definition: XLALError.c:363
int XLALIIRFilterREAL8Vector(REAL8Vector *vector, REAL8IIRFilter *filter)
double REAL8
Double precision real floating-point number (8 bytes).
#define LALMalloc(n)
Definition: LALMalloc.h:93
#define ABORT(statusptr, code, mesg)
Vector of type REAL8, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:154
REAL4 * data
Pointer to the data array.
Definition: LALDatatypes.h:150
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
void LALDIIRFilterREAL4Vector(LALStatus *stat, REAL4Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
#define IIRFILTERH_EMEM
Memory allocation error.
Definition: IIRFilter.h:133
void LALIIRFilterREAL4Vector(LALStatus *stat, REAL4Vector *vector, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
int32_t INT4
Four-byte signed integer.
void LALIIRFilterREAL8Vector(LALStatus *stat, REAL8Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
#define ASSERT(assertion, statusptr, code, mesg)
UINT4 length
Number of elements in array.
Definition: LALDatatypes.h:149
REAL4Vector * recursCoef
The recursive filter coefficients.
Definition: IIRFilter.h:156