LAL  7.2.0.1-a5e7375
IIRFilterTest.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 /** \cond DONT_DOXYGEN */
21 /** \endcond */
22 #include <complex.h>
23 #include <lal/LALStdlib.h>
24 #include <lal/LALConstants.h>
25 #include <ctype.h>
26 #include <math.h>
27 #include <string.h>
28 #include <stdlib.h>
29 #include <time.h>
30 #include <lal/AVFactories.h>
31 #include <lal/Random.h>
32 #include <lal/IIRFilter.h>
33 #include <lal/ZPGFilter.h>
34 
35 
36 /**
37  * \author Creighton, T. D.
38  * \file
39  * \ingroup IIRFilter_h
40  *
41  * \brief Tests the routines in \ref IIRFilter_h.
42  *
43  * ### Usage ###
44  *
45  * \code
46  * IIRFilterTest [-f filtertag] [-o] [-t] [-d debuglevel] [-n npts] [-r reps] [-w freq]
47  * \endcode
48  *
49  * ### Description ###
50  *
51  * This program generates a time series vector, and passes it through a
52  * third-order Butterworth low-pass filter. By default, running this
53  * program with no arguments simply passes an impulse function to the
54  * filter routines, producing no output. All filter parameters are set
55  * from <tt>\#define</tt>d constants. The following option flags are
56  * accepted:
57  * <ul>
58  * <li>[<tt>-f</tt>] Specifies which filter(s) to be used:
59  * \c filtertag is a token containing one or more character codes
60  * from \c a to \c f and/or from \c A to \c D, each
61  * corresponding to a different filter:
62  *
63  * <table>
64  * <tr><td>\c a</td><td>=</td><td><tt>LALSIIRFilter()</tt></td><td>\c A</td><td>=</td><td><tt>LALDIIRFilter()</tt></td></tr>
65  * <tr><td>\c b</td><td>=</td><td><tt>LALIIRFilterREAL4()</tt></td><td>\c B</td><td>=</td><td><tt>LALIIRFilterREAL8()</tt></td></tr>
66  * <tr><td>\c c</td><td>=</td><td><tt>LALIIRFilterREAL4Vector()</tt></td><td>\c C</td><td>=</td><td><tt>LALIIRFilterREAL8Vector()</tt></td></tr>
67  * <tr><td>\c d</td><td>=</td><td><tt>LALIIRFilterREAL4VectorR()</tt></td><td>\c D</td><td>=</td><td><tt>LALIIRFilterREAL8VectorR()</tt></td></tr>
68  * <tr><td>\c e</td><td>=</td><td><tt>LALDIIRFilterREAL4Vector()</tt></td><td></td><td></td><td></td></tr>
69  * <tr><td>\c f</td><td>=</td><td><tt>LALDIIRFilterREAL4VectorR()</tt></td><td></td><td></td><td></td></tr>
70  * </table>
71  *
72  * If not specified, <tt>-f abcd</tt> is assumed.</li>
73  *
74  * <li>[<tt>-o</tt>] Prints the input and output vectors to data files:
75  * <tt>out.0</tt> stores the initial impulse, <tt>out.</tt>\f$c\f$ the response
76  * computed using the filter with character code \f$c\f$ (above). If not
77  * specified, the routines are exercised, but no output is written.</li>
78  *
79  * <li>[<tt>-t</tt>] Causes \c IIRFilterTest to fill the time series
80  * vector with Gaussian random deviates, and prints execution times for
81  * the various filter subroutines to \c stdout. To generate useful
82  * timing data, the default size of the time vector is increased (unless
83  * explicitly set, below).</li>
84  *
85  * <li>[<tt>-d</tt>] Changes the debug level from 0 to the specified
86  * value \c debuglevel.</li>
87  *
88  * <li>[<tt>-n</tt>] Sets the size of the time vectors to \c npts.
89  * If not specified, 4096 points are used (4194304 if the <tt>-t</tt>
90  * option was also given).</li>
91  *
92  * <li>[<tt>-r</tt>] Applies each filter to the data \c reps times
93  * instead of just once.</li>
94  *
95  * <li>[<tt>-w</tt>] Sets the characteristic frequency of the filter to
96  * \c freq in the \f$w\f$-plane (described in \ref ZPGFilter.h). If
97  * not specified, <tt>-w 0.01</tt> is assumed (i.e.\ a characteristic
98  * frequency of 2\% of Nyquist).</li>
99  * </ul>
100  *
101  * ### Algorithm ###
102  *
103  * \anchor tdfilter_iirfiltertest
104  * \image html tdfilter_iirfiltertest.png "Impulse response functions in the frequency domain (computed as a windowed FFT) for various filtering routines\, run using the <tt>-r 5</tt> option."
105  *
106  * A third-order Butterworth low-pass filter is defined by the following
107  * power response function:
108  * \f[
109  * |T(w)|^2 = \frac{1}{1+(w/w_c)^6}\;,
110  * \f]
111  * where the frequency parameter \f$w=\tan(\pi f\Delta t)\f$ maps the Nyquist
112  * interval onto the entire real axis, and \f$w_c\f$ is the (transformed)
113  * characteristic frequency set by the <tt>-w</tt> option above. A stable
114  * time-domain filter has poles with \f$|z|<1\f$ in the \f$z\f$-plane
115  * representation, which means that the \f$w\f$-plane poles should have
116  * \f$\Im(w)>0\f$. We construct a transfer function using only the
117  * positive-imaginary poles of the power response function:
118  * \f[
119  * T(w) = \frac{iw_c^3}{\prod_{k=0}^2
120  * \left[w-w_c e^{(2k+1)i\pi/6}\right]}\;,
121  * \f]
122  * where we have chosen a phase coefficient \f$i\f$ in the numerator in order
123  * to get a purely real DC response. This ensures that the \f$z\f$-plane
124  * representation of the filter will have a real gain, resulting in a
125  * physically-realizable IIR filter.
126  *
127  * The poles and gain of the transfer function \f$T(w)\f$ are simply read off
128  * of the equation above, and are stored in a \c COMPLEX8ZPGFilter.
129  * This is transformed from the \f$w\f$-plane to the \f$z\f$-plane representation
130  * using <tt>LALWToZCOMPLEX8ZPGFilter()</tt>, and then used to create an
131  * IIR filter with <tt>LALCreateREAL4IIRFilter()</tt>. This in turn is
132  * used by the routines <tt>LALSIIRFilter()</tt>,
133  * <tt>LALIIRFilterREAL4()</tt>, <tt>LALIIRFilterREAL4Vector()</tt>, and
134  * <tt>LALIIRFilterREAL4VectorR()</tt> to filter a data vector containing
135  * either a unit impulse or white Gaussian noise (for more useful timing
136  * information).
137  *
138  * ### Sample output ###
139  *
140  * Running this program on a 1.3 GHz Intel machine with no optimization
141  * produced the following typical timing information:
142  * \code
143  * > IIRFilterTest -r 5 -t -f abcdefABCD
144  * Filtering 4194304 points 5 times:
145  * Elapsed time for LALSIIRFilter(): 1.39 s
146  * Elapsed time for LALDIIRFilter(): 1.79 s
147  * Elapsed time for LALIIRFilterREAL4(): 2.86 s
148  * Elapsed time for LALIIRFilterREAL8(): 3.25 s
149  * Elapsed time for LALIIRFilterREAL4Vector(): 1.52 s
150  * Elapsed time for LALIIRFilterREAL8Vector(): 2.13 s
151  * Elapsed time for LALIIRFilterREAL4VectorR(): 1.33 s
152  * Elapsed time for LALIIRFilterREAL8VectorR(): 1.96 s
153  * Elapsed time for LALDIIRFilterREAL4Vector(): 1.12 s
154  * Elapsed time for LALDIIRFilterREAL4VectorR(): 1.06 s
155  * \endcode
156  * From these results it is clear that the mixed-precision vector
157  * filtering routines are the most efficient, outperforming even the
158  * purely single-precision vector filtering routines by 20\%--30\%. This
159  * was unanticipated; by my count the main inner loop of the
160  * single-precision routines contain \f$2M+2N+1\f$ dereferences and \f$2M+2N-3\f$
161  * floating-point operations per vector element, where \f$M\f$ and \f$N\f$ are
162  * the direct and recursive filter orders, whereas the mixed-precision
163  * routines contain \f$2\max\{M,N\}+2M+2\f$ dereferences and \f$2M+2N-1\f$
164  * floating-point operations per element. However, most of the
165  * dereferences in the mixed-precision routines are to short internal
166  * arrays rather than to the larger data vector, which might cause some
167  * speedup.
168  *
169  * Running the same command with the <tt>-t</tt> flag replaced with
170  * <tt>-o</tt> generates files containing the impulse response of the
171  * filters. The frequency-domain impulse response is shown in
172  * \ref tdfilter_iirfiltertest "this figure". This shows the steady improvement in
173  * truncation error from single- to mixed- to double-precision filtering.
174  *
175  */
176 
177 /** \name Error Codes */
178 /** @{ */
179 #define IIRFILTERTESTC_ENORM 0 /**< Normal exit */
180 #define IIRFILTERTESTC_ESUB 1 /**< Subroutine failed */
181 #define IIRFILTERTESTC_EARG 2 /**< Error parsing arguments */
182 #define IIRFILTERTESTC_EBAD 3 /**< Bad argument value */
183 #define IIRFILTERTESTC_EFILE 4 /**< Could not create output file */
184 /** @} */
185 
186 /** \cond DONT_DOXYGEN */
187 #define IIRFILTERTESTC_MSGENORM "Normal exit"
188 #define IIRFILTERTESTC_MSGESUB "Subroutine failed"
189 #define IIRFILTERTESTC_MSGEARG "Error parsing arguments"
190 #define IIRFILTERTESTC_MSGEBAD "Bad argument value"
191 #define IIRFILTERTESTC_MSGEFILE "Could not create output file"
192 
193 /* Default parameters. */
194 #define NPTS 4096 /* Default length of time series */
195 #define NPTS_T 4194304 /* Length of time series for timing runs */
196 #define WC 0.01 /* Characteristic frequency in w-plane */
197 #define REPS 1 /* Number of repeated filterings */
198 #define TAG "abcd" /* Filter codes */
199 
200 /* Mathematical constant. */
201 #define SQRT3_2 0.8660254037844386467637231707529361L
202 
203 /* Usage format string. */
204 #define USAGE "Usage: %s [-d debuglevel] [-f filtertag] [-n npts]\n" \
205 "\t[-r reps] [-w freq] [-o] [-t]\n"
206 
207 /* Macros for printing errors and testing subroutines. */
208 #define ERROR( code, msg, statement ) \
209 do { \
210  if ( lalDebugLevel & LALERROR ) \
211  LALPrintError( "Error[0] %d: program %s, file %s, line %d, %s\n" \
212  " %s %s\n", (code), *argv, __FILE__, \
213  __LINE__, "$Id$", statement ? statement : \
214  "", (msg) ); \
215 } while (0)
216 
217 #define INFO( statement ) \
218 do { \
219  if ( lalDebugLevel & LALINFO ) \
220  LALPrintError( "Info[0]: program %s, file %s, line %d, %s\n" \
221  " %s\n", *argv, __FILE__, __LINE__, \
222  "$Id$", (statement) ); \
223 } while (0)
224 
225 #define SUB( func, statusptr ) \
226 do { \
227  if ( (func), (statusptr)->statusCode ) { \
228  ERROR( IIRFILTERTESTC_ESUB, IIRFILTERTESTC_MSGESUB, \
229  "Function call \"" #func "\" failed:" ); \
230  return IIRFILTERTESTC_ESUB; \
231  } \
232 } while (0)
233 
234 /* A script to assign the zpgFilter, which may be single or double
235  precision. */
236 #define ASSIGNFILTER \
237 do { \
238  zpgFilter->poles->data[0] = wc*SQRT3_2; \
239  zpgFilter->poles->data[0] += I*wc*0.5; \
240  zpgFilter->poles->data[1] = 0.0; \
241  zpgFilter->poles->data[1] += I*wc; \
242  zpgFilter->poles->data[2] = -wc*SQRT3_2; \
243  zpgFilter->poles->data[2] += I*wc*0.5; \
244  zpgFilter->gain = 0.0; \
245  zpgFilter->gain = I*wc*wc*wc; \
246 } while (0)
247 
248 /* A script to print a data vector, which may be single or double
249  precision. */
250 #define PRINTDATA( outfile ) \
251 do { \
252  FILE *fp = fopen( (outfile), "w" ); \
253  if ( !fp ) { \
254  ERROR( IIRFILTERTESTC_EFILE, IIRFILTERTESTC_MSGEFILE, 0 ); \
255  return IIRFILTERTESTC_EFILE; \
256  } \
257  i = npts; \
258  while ( i-- ) \
259  fprintf( fp, "%23.16e\n", *(data++) ); \
260 } while (0)
261 
262 
263 int
264 main(int argc, char **argv)
265 {
266  static LALStatus stat; /* Status pointer for subroutines */
267  BOOLEAN print = 0; /* Whether output will be printed */
268  BOOLEAN doTime = 0; /* Whether filters will be timed */
269  BOOLEAN points = 0; /* Whether number of points is given */
270  const CHAR *tag = TAG; /* List of filter character codes */
271  INT4 npts = NPTS; /* Number of points in time series */
272  INT4 reps = REPS; /* Number of repeated filterings */
273  INT4 arg; /* argument counter */
274  INT4 i, j; /* Index counters */
275  REAL8 wc = WC; /* Characteristic w-plane frequency */
276  REAL4Vector *sInput = NULL; /* REAL4 time series input vector */
277  REAL8Vector *dInput = NULL; /* REAL8 time series input vector */
278  REAL4Vector *sOutput = NULL; /* REAL4 time series output vector */
279  REAL8Vector *dOutput = NULL; /* REAL8 time series output vector */
280  REAL4IIRFilter *sIIRFilter = NULL; /* The REAL4 IIR filter */
281  REAL8IIRFilter *dIIRFilter = NULL; /* The REAL8 IIR filter */
282  clock_t start; /* Clock time before starting to filter */
283  clock_t stop; /* Clock time after filtering */
284 
285  /* Parse argument list. i stores the current position. */
286  arg = 1;
287  while ( arg < argc ) {
288  /* Parse debuglevel option. */
289  if ( !strcmp( argv[arg], "-d" ) ) {
290  if ( argc > arg + 1 ) {
291  arg++;
292  } else {
293  ERROR( IIRFILTERTESTC_EARG, IIRFILTERTESTC_MSGEARG, 0 );
294  LALPrintError( USAGE, *argv );
295  return IIRFILTERTESTC_EARG;
296  }
297  }
298  /* Parse filtertag option. */
299  else if ( !strcmp( argv[arg], "-f" ) ) {
300  if ( argc > arg + 1 ) {
301  arg++;
302  tag = argv[arg++];
303  } else {
304  ERROR( IIRFILTERTESTC_EARG, IIRFILTERTESTC_MSGEARG, 0 );
305  LALPrintError( USAGE, *argv );
306  return IIRFILTERTESTC_EARG;
307  }
308  }
309  /* Parse number of points. */
310  else if ( !strcmp( argv[arg], "-n" ) ) {
311  if ( argc > arg + 1 ) {
312  arg++;
313  npts = atoi( argv[arg++] );
314  points = 1;
315  } else {
316  ERROR( IIRFILTERTESTC_EARG, IIRFILTERTESTC_MSGEARG, 0 );
317  LALPrintError( USAGE, *argv );
318  return IIRFILTERTESTC_EARG;
319  }
320  }
321  /* Parse number of repetitions. */
322  else if ( !strcmp( argv[arg], "-r" ) ) {
323  if ( argc > arg + 1 ) {
324  arg++;
325  reps = strtoul( argv[arg++], NULL, 10 );
326  } else {
327  ERROR( IIRFILTERTESTC_EARG, IIRFILTERTESTC_MSGEARG, 0 );
328  LALPrintError( USAGE, *argv );
329  return IIRFILTERTESTC_EARG;
330  }
331  }
332  /* Parse characteristic frequency. */
333  else if ( !strcmp( argv[arg], "-w" ) ) {
334  if ( argc > arg + 1 ) {
335  arg++;
336  wc = atof( argv[arg++] );
337  } else {
338  ERROR( IIRFILTERTESTC_EARG, IIRFILTERTESTC_MSGEARG, 0 );
339  LALPrintError( USAGE, *argv );
340  return IIRFILTERTESTC_EARG;
341  }
342  }
343  /* Parse timing and output flags. */
344  else if ( !strcmp( argv[arg], "-t" ) ) {
345  arg++;
346  doTime = 1;
347  }
348  else if ( !strcmp( argv[arg], "-o" ) ) {
349  arg++;
350  print = 1;
351  }
352  /* Unrecognized option. */
353  else {
354  ERROR( IIRFILTERTESTC_EARG, IIRFILTERTESTC_MSGEARG, 0 );
355  LALPrintError( USAGE, *argv );
356  return IIRFILTERTESTC_EARG;
357  }
358  } /* End of argument parsing loop. */
359 
360  /* Check argument values. */
361  if ( doTime && !points )
362  npts = NPTS_T;
363  if ( npts <= 0 ) {
364  ERROR( IIRFILTERTESTC_EBAD, IIRFILTERTESTC_MSGEBAD, "npts <= 0:" );
365  LALPrintError( USAGE, *argv );
366  return IIRFILTERTESTC_EBAD;
367  }
368  if ( reps <= 0 ) {
369  ERROR( IIRFILTERTESTC_EBAD, IIRFILTERTESTC_MSGEBAD, "reps <= 0:" );
370  LALPrintError( USAGE, *argv );
371  return IIRFILTERTESTC_EBAD;
372  }
373  if ( wc <= 0.0 ) {
374  ERROR( IIRFILTERTESTC_EBAD, IIRFILTERTESTC_MSGEBAD, "freq <= 0:" );
375  LALPrintError( USAGE, *argv );
376  return IIRFILTERTESTC_EBAD;
377  }
378 
379  /* Create the time-domain filter(s). */
380  if ( strchr( tag, 'a' ) || strchr( tag, 'b' ) ||
381  strchr( tag, 'c' ) || strchr( tag, 'd' ) ) {
382  COMPLEX8ZPGFilter *zpgFilter = NULL;
383  SUB( LALCreateCOMPLEX8ZPGFilter( &stat, &zpgFilter, 0, 3 ),
384  &stat );
385  ASSIGNFILTER;
386  SUB( LALWToZCOMPLEX8ZPGFilter( &stat, zpgFilter ), &stat );
387  SUB( LALCreateREAL4IIRFilter( &stat, &sIIRFilter, zpgFilter ),
388  &stat );
389  SUB( LALDestroyCOMPLEX8ZPGFilter( &stat, &zpgFilter ), &stat );
390  }
391  if ( strchr( tag, 'A' ) || strchr( tag, 'B' ) ||
392  strchr( tag, 'C' ) || strchr( tag, 'D' ) ||
393  strchr( tag, 'e' ) || strchr( tag, 'f' ) ) {
394  COMPLEX16ZPGFilter *zpgFilter = NULL;
395  SUB( LALCreateCOMPLEX16ZPGFilter( &stat, &zpgFilter, 0, 3 ),
396  &stat );
397  ASSIGNFILTER;
398  SUB( LALWToZCOMPLEX16ZPGFilter( &stat, zpgFilter ), &stat );
399  SUB( LALCreateREAL8IIRFilter( &stat, &dIIRFilter, zpgFilter ),
400  &stat );
401  SUB( LALDestroyCOMPLEX16ZPGFilter( &stat, &zpgFilter ), &stat );
402  }
403  if ( !sIIRFilter && !dIIRFilter ) {
404  ERROR( IIRFILTERTESTC_EBAD, IIRFILTERTESTC_MSGEBAD, "filtertag:" );
405  LALPrintError( USAGE, *argv );
406  return IIRFILTERTESTC_EBAD;
407  }
408 
409  /* Create the input time series. */
410  SUB( LALSCreateVector( &stat, &sInput, npts ), &stat );
411  if ( doTime ) {
412  RandomParams *params; /* Params for random generator */
413  params = XLALCreateRandomParams( 0 );
414  INFO( "Begining generation of Gaussian random deviates.\n" );
415  XLALNormalDeviates( sInput, params );
416  INFO( "Finished generating random deviates.\n" );
417  XLALDestroyRandomParams( params );
418  } else {
419  memset( sInput->data, 0, npts*sizeof(REAL4) );
420  sInput->data[npts/2] = 1.0;
421  }
422  if ( print ) {
423  REAL4 *data = sInput->data;
424  PRINTDATA( "out.0" );
425  }
426 
427  /* Create other time series. */
428  if ( strchr( tag, 'a' ) || strchr( tag, 'b' ) ||
429  strchr( tag, 'c' ) || strchr( tag, 'd' ) ||
430  strchr( tag, 'e' ) || strchr( tag, 'f' ) ) {
431  SUB( LALSCreateVector( &stat, &sOutput, npts ), &stat );
432  }
433  if ( strchr( tag, 'A' ) || strchr( tag, 'B' ) ||
434  strchr( tag, 'C' ) || strchr( tag, 'D' ) ) {
435  SUB( LALDCreateVector( &stat, &dInput, npts ), &stat );
436  SUB( LALDCreateVector( &stat, &dOutput, npts ), &stat );
437  for ( i = 0; i < npts; i++ )
438  dInput->data[i] = (REAL8)( sInput->data[i] );
439  }
440 
441  /* Filter the time series. */
442  if ( doTime ) {
443  if ( reps == 1 )
444  fprintf( stdout, "Filtering %i points:\n", npts );
445  else
446  fprintf( stdout, "Filtering %i points %i times:\n", npts,
447  reps );
448  }
449 
450  /* Using LALSIIRFilter(): */
451  if ( strchr( tag, 'a' ) ) {
452  memcpy( sOutput->data, sInput->data, npts*sizeof(REAL4) );
453  j = reps;
454  start = clock();
455  while ( j-- ) {
456  REAL4 *data = sOutput->data;
457  i = npts;
458  while ( i-- ) {
459  *data = LALSIIRFilter( *data, sIIRFilter );
460  data++;
461  }
462  }
463  stop = clock();
464  if ( doTime )
465  fprintf( stdout, "Elapsed time for LALSIIRFilter(): %.2f"
466  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
467  if ( print ) {
468  REAL4 *data = sOutput->data;
469  PRINTDATA( "out.a" );
470  }
471  }
472 
473  /* Using LALDIIRFilter(): */
474  if ( strchr( tag, 'A' ) ) {
475  memcpy( dOutput->data, dInput->data, npts*sizeof(REAL8) );
476  j = reps;
477  start = clock();
478  while ( j-- ) {
479  REAL8 *data = dOutput->data;
480  i = npts;
481  while ( i-- ) {
482  *data = LALDIIRFilter( *data, dIIRFilter );
483  data++;
484  }
485  }
486  stop = clock();
487  if ( doTime )
488  fprintf( stdout, "Elapsed time for LALDIIRFilter(): %.2f"
489  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
490  if ( print ) {
491  REAL8 *data = dOutput->data;
492  PRINTDATA( "out.A" );
493  }
494  }
495 
496  /* Using LALIIRFilterREAL4(). */
497  if ( strchr( tag, 'b' ) ) {
498  memset( sIIRFilter->history->data, 0,
499  sIIRFilter->history->length*sizeof(REAL4) );
500  memcpy( sOutput->data, sInput->data, npts*sizeof(REAL4) );
501  j = reps;
502  start = clock();
503  while ( j-- ) {
504  REAL4 *data = sOutput->data;
505  i = npts;
506  while ( i-- ) {
507  SUB( LALIIRFilterREAL4( &stat, data, *data, sIIRFilter ),
508  &stat );
509  data++;
510  }
511  }
512  stop = clock();
513  if ( doTime )
514  fprintf( stdout, "Elapsed time for LALIIRFilterREAL4(): %.2f"
515  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
516  if( print ) {
517  REAL4 *data = sOutput->data;
518  PRINTDATA( "out.b" );
519  }
520  }
521 
522  /* Using LALIIRFilterREAL8(). */
523  if ( strchr( tag, 'B' ) ) {
524  memset( dIIRFilter->history->data, 0,
525  dIIRFilter->history->length*sizeof(REAL8) );
526  memcpy( dOutput->data, dInput->data, npts*sizeof(REAL8) );
527  j = reps;
528  start = clock();
529  while ( j-- ) {
530  REAL8 *data = dOutput->data;
531  i = npts;
532  while ( i-- ) {
533  SUB( LALIIRFilterREAL8( &stat, data, *data, dIIRFilter ),
534  &stat );
535  data++;
536  }
537  }
538  stop = clock();
539  if ( doTime )
540  fprintf( stdout, "Elapsed time for LALIIRFilterREAL8(): %.2f"
541  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
542  if( print ) {
543  REAL8 *data = dOutput->data;
544  PRINTDATA( "out.B" );
545  }
546  }
547 
548  /* Using LALIIRFilterREAL4Vector(). */
549  if ( strchr( tag, 'c' ) ) {
550  memset( sIIRFilter->history->data, 0,
551  sIIRFilter->history->length*sizeof(REAL4) );
552  memcpy( sOutput->data, sInput->data, npts*sizeof(REAL4) );
553  j = reps;
554  start = clock();
555  while ( j-- )
556  SUB( LALIIRFilterREAL4Vector( &stat, sOutput, sIIRFilter ),
557  &stat );
558  stop=clock();
559  if( doTime )
560  fprintf( stdout, "Elapsed time for LALIIRFilterREAL4Vector(): %.2f"
561  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
562  if ( print ) {
563  REAL4 *data = sOutput->data;
564  PRINTDATA( "out.c" );
565  }
566  }
567 
568  /* Using LALIIRFilterREAL8Vector(). */
569  if ( strchr( tag, 'C' ) ) {
570  memset( dIIRFilter->history->data, 0,
571  dIIRFilter->history->length*sizeof(REAL8) );
572  memcpy( dOutput->data, dInput->data, npts*sizeof(REAL8) );
573  j = reps;
574  start = clock();
575  while ( j-- )
576  SUB( LALIIRFilterREAL8Vector( &stat, dOutput, dIIRFilter ),
577  &stat );
578  stop=clock();
579  if( doTime )
580  fprintf( stdout, "Elapsed time for LALIIRFilterREAL8Vector(): %.2f"
581  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
582  if ( print ) {
583  REAL8 *data = dOutput->data;
584  PRINTDATA( "out.C" );
585  }
586  }
587 
588  /* Using LALIIRFilterREAL4VectorR(). */
589  if ( strchr( tag, 'd' ) ) {
590  memcpy( sOutput->data, sInput->data, npts*sizeof(REAL4) );
591  j = reps;
592  start = clock();
593  while ( j-- )
594  SUB( LALIIRFilterREAL4VectorR( &stat, sOutput, sIIRFilter ),
595  &stat );
596  stop=clock();
597  if ( doTime )
598  fprintf( stdout, "Elapsed time for LALIIRFilterREAL4VectorR(): %.2f"
599  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
600  if( print ) {
601  REAL4 *data = sOutput->data;
602  PRINTDATA( "out.d" );
603  }
604  }
605 
606  /* Using LALIIRFilterREAL8VectorR(). */
607  if ( strchr( tag, 'D' ) ) {
608  memcpy( dOutput->data, dInput->data, npts*sizeof(REAL8) );
609  j = reps;
610  start = clock();
611  while ( j-- )
612  SUB( LALIIRFilterREAL8VectorR( &stat, dOutput, dIIRFilter ),
613  &stat );
614  stop=clock();
615  if ( doTime )
616  fprintf( stdout, "Elapsed time for LALIIRFilterREAL8VectorR(): %.2f"
617  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
618  if( print ) {
619  REAL8 *data = dOutput->data;
620  PRINTDATA( "out.D" );
621  }
622  }
623 
624  /* Using LALDIIRFilterREAL4Vector(). */
625  if ( strchr( tag, 'e' ) ) {
626  memset( dIIRFilter->history->data, 0,
627  dIIRFilter->history->length*sizeof(REAL8) );
628  memcpy( sOutput->data, sInput->data, npts*sizeof(REAL4) );
629  j = reps;
630  start = clock();
631  while ( j-- )
632  SUB( LALDIIRFilterREAL4Vector( &stat, sOutput, dIIRFilter ),
633  &stat );
634  stop=clock();
635  if( doTime )
636  fprintf( stdout, "Elapsed time for LALDIIRFilterREAL4Vector(): %.2f"
637  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
638  if ( print ) {
639  REAL4 *data = sOutput->data;
640  PRINTDATA( "out.e" );
641  }
642  }
643 
644  /* Using LALDIIRFilterREAL4VectorR(). */
645  if ( strchr( tag, 'f' ) ) {
646  memcpy( sOutput->data, sInput->data, npts*sizeof(REAL4) );
647  j = reps;
648  start = clock();
649  while ( j-- )
650  SUB( LALDIIRFilterREAL4VectorR( &stat, sOutput, dIIRFilter ),
651  &stat );
652  stop=clock();
653  if( doTime )
654  fprintf( stdout, "Elapsed time for LALDIIRFilterREAL4VectorR(): %.2f"
655  " s\n", (double)( stop - start )/CLOCKS_PER_SEC );
656  if ( print ) {
657  REAL4 *data = sOutput->data;
658  PRINTDATA( "out.f" );
659  }
660  }
661 
662  /* Free memory and exit. */
663  SUB( LALSDestroyVector( &stat, &sInput ), &stat );
664  if ( sOutput )
665  SUB( LALSDestroyVector( &stat, &sOutput ), &stat );
666  if ( dInput )
667  SUB( LALDDestroyVector( &stat, &dInput ), &stat );
668  if ( dOutput )
669  SUB( LALDDestroyVector( &stat, &dOutput ), &stat );
670  if ( sIIRFilter )
671  SUB( LALDestroyREAL4IIRFilter( &stat, &sIIRFilter ), &stat );
672  if ( dIIRFilter )
673  SUB( LALDestroyREAL8IIRFilter( &stat, &dIIRFilter ), &stat );
675  INFO( IIRFILTERTESTC_MSGENORM );
676  return IIRFILTERTESTC_ENORM;
677 }
678 /** \endcond */
void LALIIRFilterREAL4VectorR(LALStatus *status, REAL4Vector *vector, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
void LALDCreateVector(LALStatus *, REAL8Vector **, UINT4)
This structure stores the direct and recursive REAL8 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:168
RandomParams * XLALCreateRandomParams(INT4 seed)
Definition: Random.c:102
void LALIIRFilterREAL8(LALStatus *stat, REAL8 *output, REAL8 input, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
Definition: IIRFilter.c:112
void LALDestroyCOMPLEX8ZPGFilter(LALStatus *stat, COMPLEX8ZPGFilter **input)
void LALSCreateVector(LALStatus *, REAL4Vector **, UINT4)
void LALCreateCOMPLEX16ZPGFilter(LALStatus *stat, COMPLEX16ZPGFilter **output, INT4 numZeros, INT4 numPoles)
#define SUB(func, statusptr)
REAL8 * data
Pointer to the data array.
Definition: LALDatatypes.h:159
#define LALDIIRFilter(x, f)
Definition: IIRFilter.h:216
void LALIIRFilterREAL8VectorR(LALStatus *status, REAL8Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
char CHAR
One-byte signed integer, see Headers LAL(Atomic)Datatypes.h for more details.
This structure stores the direct and recursive REAL4 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:152
void LALCheckMemoryLeaks(void)
Definition: LALMalloc.c:784
See DATATYPE-ZPGFilter types for details.
Definition: LALDatatypes.h:930
REAL8Vector * history
The previous values of w.
Definition: IIRFilter.h:173
void LALDIIRFilterREAL4VectorR(LALStatus *status, REAL4Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
See DATATYPE-ZPGFilter types for details.
Definition: LALDatatypes.h:921
REAL4Vector * history
The previous values of w.
Definition: IIRFilter.h:157
#define fprintf
float REAL4
Single precision real floating-point number (4 bytes).
#define IIRFILTERTESTC_ENORM
Normal exit.
Vector of type REAL4, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:145
#define IIRFILTERTESTC_EARG
Error parsing arguments.
void LALWToZCOMPLEX8ZPGFilter(LALStatus *stat, COMPLEX8ZPGFilter *filter)
Deprecated.
int XLALNormalDeviates(REAL4Vector *deviates, RandomParams *params)
Definition: Random.c:171
void LALDestroyREAL4IIRFilter(LALStatus *stat, REAL4IIRFilter **input)
Deprecated.
void LALSDestroyVector(LALStatus *, REAL4Vector **)
void LALCreateREAL8IIRFilter(LALStatus *stat, REAL8IIRFilter **output, COMPLEX16ZPGFilter *input)
Deprecated.
void LALWToZCOMPLEX16ZPGFilter(LALStatus *stat, COMPLEX16ZPGFilter *filter)
Deprecated.
double REAL8
Double precision real floating-point number (8 bytes).
void LALCreateCOMPLEX8ZPGFilter(LALStatus *stat, COMPLEX8ZPGFilter **output, INT4 numZeros, INT4 numPoles)
#define IIRFILTERTESTC_EBAD
Bad argument value.
Vector of type REAL8, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:154
void LALIIRFilterREAL4(LALStatus *stat, REAL4 *output, REAL4 input, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
Definition: IIRFilter.c:55
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 LALDestroyREAL8IIRFilter(LALStatus *stat, REAL8IIRFilter **input)
Deprecated.
void LALDIIRFilterREAL4Vector(LALStatus *status, REAL4Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
#define INFO(statement)
void LALDestroyCOMPLEX16ZPGFilter(LALStatus *stat, COMPLEX16ZPGFilter **input)
unsigned char BOOLEAN
Boolean logical type, see Headers LAL(Atomic)Datatypes.h for more details.
UINT4 length
Number of elements in array.
Definition: LALDatatypes.h:158
void LALIIRFilterREAL4Vector(LALStatus *status, REAL4Vector *vector, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
#define ERROR(code, msg, statement)
int main(int argc, char *argv[])
Definition: cache.c:25
int32_t INT4
Four-byte signed integer.
#define USAGE
Definition: GzipTest.c:27
void LALDDestroyVector(LALStatus *, REAL8Vector **)
void LALIIRFilterREAL8Vector(LALStatus *status, REAL8Vector *vector, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
REAL4 LALSIIRFilter(REAL4 x, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
Definition: IIRFilter.c:143
void LALCreateREAL4IIRFilter(LALStatus *stat, REAL4IIRFilter **output, COMPLEX8ZPGFilter *input)
Deprecated.
UINT4 length
Number of elements in array.
Definition: LALDatatypes.h:149
int LALPrintError(const char *fmt,...)
Definition: LALError.c:46
void XLALDestroyRandomParams(RandomParams *params)
Definition: Random.c:140
This structure contains the parameters necessary for generating the current sequence of random number...
Definition: Random.h:84