Loading [MathJax]/extensions/TeX/AMSsymbols.js
LAL 7.7.0.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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 ) \
209do { \
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 ) \
218do { \
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 ) \
226do { \
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 \
237do { \
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 ) \
251do { \
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
263int
264main(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 );
677}
678/** \endcond */
#define USAGE
Definition: GzipTest.c:27
#define LALDIIRFilter(x, f)
Definition: IIRFilter.h:216
#define IIRFILTERTESTC_EBAD
Bad argument value.
#define IIRFILTERTESTC_ENORM
Normal exit.
#define IIRFILTERTESTC_EARG
Error parsing arguments.
void LALCheckMemoryLeaks(void)
Definition: LALMalloc.c:784
#define SUB(func, statusptr)
#define ERROR(code, msg, statement)
#define INFO(statement)
#define fprintf
int main(int argc, char *argv[])
Definition: cache.c:25
void LALWToZCOMPLEX8ZPGFilter(LALStatus *stat, COMPLEX8ZPGFilter *filter)
Deprecated.
void LALWToZCOMPLEX16ZPGFilter(LALStatus *stat, COMPLEX16ZPGFilter *filter)
Deprecated.
void LALCreateREAL4IIRFilter(LALStatus *stat, REAL4IIRFilter **output, COMPLEX8ZPGFilter *input)
Deprecated.
void LALCreateREAL8IIRFilter(LALStatus *stat, REAL8IIRFilter **output, COMPLEX16ZPGFilter *input)
Deprecated.
void LALCreateCOMPLEX8ZPGFilter(LALStatus *stat, COMPLEX8ZPGFilter **output, INT4 numZeros, INT4 numPoles)
void LALCreateCOMPLEX16ZPGFilter(LALStatus *stat, COMPLEX16ZPGFilter **output, INT4 numZeros, INT4 numPoles)
void LALDestroyREAL4IIRFilter(LALStatus *stat, REAL4IIRFilter **input)
Deprecated.
void LALDestroyREAL8IIRFilter(LALStatus *stat, REAL8IIRFilter **input)
Deprecated.
void LALDestroyCOMPLEX16ZPGFilter(LALStatus *stat, COMPLEX16ZPGFilter **input)
void LALDestroyCOMPLEX8ZPGFilter(LALStatus *stat, COMPLEX8ZPGFilter **input)
void LALIIRFilterREAL8(LALStatus *stat, REAL8 *output, REAL8 input, REAL8IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
Definition: IIRFilter.c:112
REAL4 LALSIIRFilter(REAL4 x, REAL4IIRFilter *filter)
WARNING: THIS FUNCTION IS OBSOLETE.
Definition: IIRFilter.c:143
void LALIIRFilterREAL4(LALStatus *stat, 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.
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).
char CHAR
One-byte signed integer, see Headers LAL(Atomic)Datatypes.h for more details.
int32_t INT4
Four-byte signed integer.
float REAL4
Single precision real floating-point number (4 bytes).
int LALPrintError(const char *fmt,...)
Definition: LALError.c:46
RandomParams * XLALCreateRandomParams(INT4 seed)
Definition: Random.c:102
void XLALDestroyRandomParams(RandomParams *params)
Definition: Random.c:140
int XLALNormalDeviates(REAL4Vector *deviates, RandomParams *params)
Definition: Random.c:171
void LALDCreateVector(LALStatus *, REAL8Vector **, UINT4)
void LALDDestroyVector(LALStatus *, REAL8Vector **)
void LALSDestroyVector(LALStatus *, REAL4Vector **)
void LALSCreateVector(LALStatus *, REAL4Vector **, UINT4)
See DATATYPE-ZPGFilter types for details.
Definition: LALDatatypes.h:930
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
Vector of type REAL4, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:145
REAL4 * data
Pointer to the data array.
Definition: LALDatatypes.h:150
UINT4 length
Number of elements in array.
Definition: LALDatatypes.h:149
This structure stores the direct and recursive REAL8 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:168
REAL8Vector * history
The previous values of w.
Definition: IIRFilter.h:173
Vector of type REAL8, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:154
REAL8 * data
Pointer to the data array.
Definition: LALDatatypes.h:159
UINT4 length
Number of elements in array.
Definition: LALDatatypes.h:158
This structure contains the parameters necessary for generating the current sequence of random number...
Definition: Random.h:86