LAL  7.1.7.1-2d066e5
CreateIIRFilter.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 
27 /**
28  * \addtogroup CreateIIRFilter_c
29  * \author Creighton, T. D.
30  *
31  * \brief Creates IIR filter objects.
32  *
33  * ### Description ###
34  *
35  * These functions create an object <tt>**output</tt> of type
36  * <tt><datatype>IIRFilter</tt>, where <tt><datatype></tt> is \c REAL4 or
37  * \c REAL8. The filter coefficients are computed from the zeroes,
38  * poles, and gain of an input object <tt>*input</tt> of type
39  * \c COMPLEX8ZPGFilter or \c COMPLEX16ZPGFilter, respectively;
40  * the sampling time interval is taken directly from
41  * <tt>input->deltaT</tt>. The ZPG filter should express the factored
42  * transfer function in the \f$z=\exp(2\pi if)\f$ plane. Initially the
43  * output handle must be a valid handle (\c output\f$\neq\f$\c NULL)
44  * but should not point to an existing object
45  * (<tt>*output</tt>=\c NULL)
46  *
47  * ### Algorithm ###
48  *
49  * An IIR filter is a real time-domain filter, which imposes certain
50  * constraints on the zeros, poles, and gain of the filter transfer
51  * function. The function <tt>Create<datatype>IIRFilter()</tt> deals with
52  * the constraints either by aborting if they are not met, or by
53  * adjusting the filter response so that they are met. In the latter
54  * case, warning messages will be issued if the external parameter
55  * \c lalDebugLevel is set to allow such messages. The specific
56  * constraints, and how they are dealt with, are as follows:
57  *
58  * First, the filter must be \e causal; that is, the output at any
59  * time can be generated entirely from the input at previous times. In
60  * practice this means that the number of (finite) poles in the \f$z\f$ plane
61  * must equal or exceed the number of (finite) zeros. If this is not the
62  * case, <tt>Create<datatype>IIRFilter()</tt> will add additional poles at
63  * \f$z=0\f$. Effectively this is just adding a delay to the filter response
64  * in order to make it causal.
65  *
66  * Second, the filter should be \e stable, which means that all poles
67  * should be located on or within the circle \f$|z|=1\f$. This is not
68  * enforced by <tt>Create<datatype>IIRFilter()</tt>, which can be used to
69  * make unstable filters; however, warnings will be issued. (In some
70  * sense the first condition is a special case of this one, since a
71  * transfer function with more zeros than poles actually has
72  * corresponding poles at infinity.)
73  *
74  * Third, the filter must be <em>physically realizable</em>; that is, the
75  * transfer function should expand to a rational function of \f$z\f$ with
76  * real coefficients. Necessary and sufficient conditions for this are
77  * that the gain be real, and that all zeros and poles either be real or
78  * come in complex conjugate pairs. The routine
79  * <tt>Create<datatype>IIRFilter()</tt> enforces this by using only the
80  * real part of the gain, and only the real or positive-imaginary zeros
81  * and poles; it assumes that the latter are paired with
82  * negative-imaginary conjugates. The routine will abort if this
83  * assumption results in a change in the given number of zeros or poles.
84  * If \c lalDebugLevel is set to allow warnings, the routine will
85  * actually check to see that each pair of nonreal poles or zeros are in
86  * fact complex conjugates, and will issue a warning if an unmatched pair
87  * is detected; however, the algorithm will then simply proceed as if the
88  * negative-imaginary points were relocated to the "correct" positions.
89  *
90  * The code associated with the warning messages is potentially rather
91  * cumbersome for production algorithms; therefore, the value of
92  * \c lalDebugLevel is tested before performing any other tests
93  * associated with warning messages. Furthermore, this code block is
94  * surrounded with compiler directives to exclude the code entirely if
95  * the module is compiled with the \c NDEBUG flag set.
96  *
97  */
98 /** @{ */
99 
100 /** \see See \ref CreateIIRFilter_c for documentation */
102 {
103  REAL4IIRFilter *output;
104  INT4 i; /* Index counter for zeros and poles. */
105  INT4 numZeros; /* The number of zeros. */
106  INT4 numPoles; /* The number of poles. */
107  INT4 numDirect; /* The number of direct filter coefficients. */
108  INT4 numRecurs; /* The number of recursive filter coefficients. */
109  INT4 num; /* An extra counter for error checking. */
110  COMPLEX8 *zeros; /* The zeros of the transfer function. */
111  COMPLEX8 *poles; /* The poles of the transfer function. */
112  REAL4 *direct; /* The direct filter coefficients. */
113  REAL4 *recurs; /* The recursive filter coefficients. */
114  REAL4 *history; /* The filter history. */
115 
116  /* Make sure all the input structures have been initialized. */
117  if ( ! input )
119  if ( ! input->zeros || ! input->poles
120  || ! input->zeros->data || ! input->poles->data )
122 
123  numZeros=input->zeros->length;
124  numPoles=input->poles->length;
125  zeros=input->zeros->data;
126  poles=input->poles->data;
127 
128  /* Check that zeros are appropriately paired. Also, keep track of
129  the number of zeros at z=0, since these will reduce the number of
130  direct coefficients required. */
131  numDirect=1;
132  for(i=0,num=0;i<numZeros;i++)
133  if(cimagf(zeros[i])==0.0){
134  num+=1;
135  if(crealf(zeros[i])==0.0)
136  numDirect-=1;
137  }
138  else if(cimagf(zeros[i])>0.0){
139  num+=2;
140 
141 #ifndef NDEBUG
143  /* Check to see that another zero is an actual conjugate.
144  This is not a foolproof test, as it can be fooled by
145  multiple zeros at the same location. */
146  INT4 j=0;
147  INT4 k=0;
148  REAL4 x = crealf(zeros[i]) - crealf(zeros[0]);
149  REAL4 y = cimagf(zeros[i]) + cimagf(zeros[0]);
150  REAL4 sep = x*x + y*y;
151  for(j=1;j<numZeros;j++){
152  x=crealf(zeros[i])-crealf(zeros[j]);
153  y=cimagf(zeros[i])+cimagf(zeros[j]);
154  if(sep>x*x+y*y){
155  sep=x*x+y*y;
156  k=j;
157  }
158  }
159  if(sep>LAL_REAL4_EPS){
160  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
161  XLALPrintWarning("Complex zero has no conjugate pair\n");
162  XLALPrintWarning("\tUnmatched zero z_%i = %.8e + i*%.8e\n",i,
163  crealf(zeros[i]),cimagf(zeros[i]));
164  XLALPrintWarning("\tNearest pair z_%i = %.8e + i*%.8e\n",k,
165  crealf(zeros[k]),cimagf(zeros[k]));
166  }
167  }
168 #endif
169  }
170 
171  if ( num != numZeros )
173 
174  /* Check that poles are appropriately paired. Also, keep track of
175  the number of poles at z=0, since these will reduce the number of
176  recursive coefficients required. */
177  numRecurs=1;
178  for(i=0,num=0;i<numPoles;i++)
179  if(cimagf(poles[i])==0.0){
180  num+=1;
181  if(crealf(poles[i])==0.0)
182  numRecurs-=1;
183  }
184  else if(cimagf(poles[i])>0.0){
185  num+=2;
186 
187 #ifndef NDEBUG
189  /* Check to see that another pole is an actual conjugate.
190  This is not a foolproof test, as it can be fooled by
191  multiple poles at the same location. */
192  INT4 j=0;
193  INT4 k=0;
194  REAL4 x = crealf(poles[i]) - crealf(poles[0]);
195  REAL4 y = cimagf(poles[i]) + cimagf(poles[0]);
196  REAL4 sep = x*x + y*y;
197  for(j=1;j<numPoles;j++){
198  x=crealf(poles[i])-crealf(poles[j]);
199  y=cimagf(poles[i])+cimagf(poles[j]);
200  if(sep>x*x+y*y){
201  sep=x*x+y*y;
202  k=j;
203  }
204  }
205  if(sep>LAL_REAL4_EPS){
206  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
207  XLALPrintWarning("Complex pole has no conjugate pair\n");
208  XLALPrintWarning("\tUnmatched pole p_%i = %.8e + i*%.8e\n",i,
209  crealf(poles[i]),cimagf(poles[i]));
210  XLALPrintWarning("\tNearest pair p_%i = %.8e + i*%.8e\n",k,
211  crealf(poles[k]),cimagf(poles[k]));
212  }
213  }
214 #endif
215  }
216 
217  if ( num != numPoles )
219 
220 #ifndef NDEBUG
222  /* Issue a warning if the gain is nonreal. */
223  if(fabs(cimag(input->gain))>fabs(LAL_REAL4_EPS*creal(input->gain))){
224  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
225  XLALPrintWarning("Gain is non-real\n");
226  XLALPrintWarning("\tg = %.8e + i*%.8e\n", crealf(input->gain), cimagf(input->gain));
227  }
228  /* Issue a warning if there are any ``removeable'' poles. */
229  for(i=0;i<numPoles;i++){
230  INT4 j=0;
231  for(;j<numZeros;j++)
232  if((crealf(poles[i])==crealf(zeros[j]))&&(cimagf(poles[i])==cimagf(zeros[j]))){
233  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
234  XLALPrintWarning("Removeable pole\n");
235  XLALPrintWarning("\tp_%i = z_%i = %.8e + i*%.8e\n",i,j,
236  crealf(poles[i]),cimagf(poles[i]));
237  }
238  }
239  /* Issue a warning if extra factors of 1/z will be applied. */
240  if(numPoles<numZeros){
241  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
242  XLALPrintWarning("Filter has more zeros than poles\n");
243  XLALPrintWarning("\t%i poles added at complex origin\n",
244  numZeros-numPoles);
245  }
246  /* Issue a warning if any poles are outside |z|=1. */
247  for(i=0;i<numPoles;i++){
248  REAL4 zAbs=crealf(poles[i])*crealf(poles[i])+cimagf(poles[i])*cimagf(poles[i]);
249  if(zAbs>1.0){
250  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
251  XLALPrintWarning("Filter has pole outside of unit circle\n");
252  XLALPrintWarning("\tp_%i = %.8e + i*%.8e, |p_%i| = %.8e\n",i,
253  crealf(poles[i]),cimagf(poles[i]),i,zAbs);
254  }
255  }
256  }
257 #endif
258 
259  /* Everything seems okay, so initialize the filter. */
260  output=LALCalloc(1,sizeof(*output));
261  if ( ! output )
263  num = (numPoles>=numZeros) ? numZeros : numPoles;
264  numDirect+=num;
265  numRecurs+=num;
266 
267  output->deltaT=input->deltaT;
268  output->directCoef = XLALCreateREAL4Vector( numDirect );
269  output->recursCoef = XLALCreateREAL4Vector( numRecurs );
270  if ( ! output->directCoef || ! output->recursCoef )
271  {
272  XLALDestroyREAL4IIRFilter( output );
274  }
275  direct=output->directCoef->data;
276  recurs=output->recursCoef->data;
277 
278  /* Expand the denominator as a polynomial in z. */
279  *recurs=-1.0;
280  for(i=1;i<numRecurs;i++)
281  recurs[i]=0.0;
282  for(i=0;i<numPoles;i++){
283  INT4 j=numRecurs-1;
284  REAL4 x=crealf(poles[i]);
285  REAL4 y=cimagf(poles[i]);
286  if(y==0.0)
287  for(;j>0;j--)
288  recurs[j]-=x*recurs[j-1];
289  else if(y>0.0){
290  for(;j>1;j--)
291  recurs[j]-=2.0*x*recurs[j-1]-(x*x+y*y)*recurs[j-2];
292  recurs[j]-=2.0*x*recurs[j-1];
293  }
294  }
295 
296  /* Expand the numerator as a polynomial in z. */
297  for(i=0;i<numDirect;i++)
298  direct[i]=0.0;
299  direct[num-numZeros]=crealf(input->gain);
300  for(i=0;i<numZeros;i++){
301  INT4 j=numDirect-1;
302  REAL4 x=crealf(zeros[i]);
303  REAL4 y=cimagf(zeros[i]);
304  if(y==0.0)
305  for(;j>num-numZeros;j--)
306  direct[j]-=x*direct[j-1];
307  else if(y>0.0){
308  for(;j>num-numZeros+1;j--)
309  direct[j]-=2.0*x*direct[j-1]-(x*x+y*y)*direct[j-2];
310  direct[j]-=2.0*x*direct[j-1];
311  }
312  }
313 
314  /* Initialize filter history. */
315  if(numDirect>=numRecurs)
316  num=numDirect-1;
317  else
318  num=numRecurs-1;
319  output->history = XLALCreateREAL4Vector( numRecurs );
320  if ( ! output->history )
321  {
322  XLALDestroyREAL4IIRFilter( output );
324  }
325  history=output->history->data;
326  for(i=0;i<num;i++)
327  history[i]=0.0;
328 
329  /* Normal exit */
330  return output;
331 }
332 
333 /** \see See \ref CreateIIRFilter_c for documentation */
335 {
336  REAL8IIRFilter *output;
337  INT4 i; /* Index counter for zeros and poles. */
338  INT4 numZeros; /* The number of zeros. */
339  INT4 numPoles; /* The number of poles. */
340  INT4 numDirect; /* The number of direct filter coefficients. */
341  INT4 numRecurs; /* The number of recursive filter coefficients. */
342  INT4 num; /* An extra counter for error checking. */
343  COMPLEX16 *zeros; /* The zeros of the transfer function. */
344  COMPLEX16 *poles; /* The poles of the transfer function. */
345  REAL8 *direct; /* The direct filter coefficients. */
346  REAL8 *recurs; /* The recursive filter coefficients. */
347  REAL8 *history; /* The filter history. */
348 
349  /* Make sure all the input structures have been initialized. */
350  if ( ! input )
352  if ( ! input->zeros || ! input->poles
353  || ! input->zeros->data || ! input->poles->data )
355 
356  numZeros=input->zeros->length;
357  numPoles=input->poles->length;
358  zeros=input->zeros->data;
359  poles=input->poles->data;
360 
361  /* Check that zeros are appropriately paired. Also, keep track of
362  the number of zeros at z=0, since these will reduce the number of
363  direct coefficients required. */
364  numDirect=1;
365  for(i=0,num=0;i<numZeros;i++)
366  if(cimag(zeros[i])==0.0){
367  num+=1;
368  if(creal(zeros[i])==0.0)
369  numDirect-=1;
370  }
371  else if(cimag(zeros[i])>0.0){
372  num+=2;
373 
374 #ifndef NDEBUG
376  /* Check to see that another zero is an actual conjugate.
377  This is not a foolproof test, as it can be fooled by
378  multiple zeros at the same location. */
379  INT4 j=0;
380  INT4 k=0;
381  REAL8 x = creal(zeros[i]) - creal(zeros[0]);
382  REAL8 y = cimag(zeros[i]) + cimag(zeros[0]);
383  REAL8 sep = x*x + y*y;
384  for(j=1;j<numZeros;j++){
385  x=creal(zeros[i])-creal(zeros[j]);
386  y=cimag(zeros[i])+cimag(zeros[j]);
387  if(sep>x*x+y*y){
388  sep=x*x+y*y;
389  k=j;
390  }
391  }
392  if(sep>LAL_REAL8_EPS){
393  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
394  XLALPrintWarning("Complex zero has no conjugate pair\n");
395  XLALPrintWarning("\tUnmatched zero z_%i = %.8e + i*%.8e\n",i,
396  creal(zeros[i]),cimag(zeros[i]));
397  XLALPrintWarning("\tNearest pair z_%i = %.8e + i*%.8e\n",k,
398  creal(zeros[k]),cimag(zeros[k]));
399  }
400  }
401 #endif
402  }
403 
404  if ( num != numZeros )
406 
407  /* Check that poles are appropriately paired. Also, keep track of
408  the number of poles at z=0, since these will reduce the number of
409  recursive coefficients required. */
410  numRecurs=1;
411  for(i=0,num=0;i<numPoles;i++)
412  if(cimag(poles[i])==0.0){
413  num+=1;
414  if(creal(poles[i])==0.0)
415  numRecurs-=1;
416  }
417  else if(cimag(poles[i])>0.0){
418  num+=2;
419 
420 #ifndef NDEBUG
422  /* Check to see that another pole is an actual conjugate.
423  This is not a foolproof test, as it can be fooled by
424  multiple poles at the same location. */
425  INT4 j=0;
426  INT4 k=0;
427  REAL8 x = creal(poles[i]) - creal(poles[0]);
428  REAL8 y = cimag(poles[i]) + cimag(poles[0]);
429  REAL8 sep = x*x + y*y;
430  for(j=1;j<numPoles;j++){
431  x=creal(poles[i])-creal(poles[j]);
432  y=cimag(poles[i])+cimag(poles[j]);
433  if(sep>x*x+y*y){
434  sep=x*x+y*y;
435  k=j;
436  }
437  }
438  if(sep>LAL_REAL8_EPS){
439  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
440  XLALPrintWarning("Complex pole has no conjugate pair\n");
441  XLALPrintWarning("\tUnmatched pole p_%i = %.8e + i*%.8e\n",i,
442  creal(poles[i]),cimag(poles[i]));
443  XLALPrintWarning("\tNearest pair p_%i = %.8e + i*%.8e\n",k,
444  creal(poles[k]),cimag(poles[k]));
445  }
446  }
447 #endif
448  }
449 
450  if ( num != numPoles )
452 
453 #ifndef NDEBUG
455  /* Issue a warning if the gain is nonreal. */
456  if(fabs(cimag(input->gain))>fabs(LAL_REAL8_EPS*creal(input->gain))){
457  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
458  XLALPrintWarning("Gain is non-real\n");
459  XLALPrintWarning("\tg = %.8e + i*%.8e\n", creal(input->gain), cimag(input->gain));
460  }
461  /* Issue a warning if there are any ``removeable'' poles. */
462  for(i=0;i<numPoles;i++){
463  INT4 j=0;
464  for(;j<numZeros;j++)
465  if((creal(poles[i])==creal(zeros[j]))&&(cimag(poles[i])==cimag(zeros[j]))){
466  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
467  XLALPrintWarning("Removeable pole\n");
468  XLALPrintWarning("\tp_%i = z_%i = %.8e + i*%.8e\n",i,j,
469  creal(poles[i]),cimag(poles[i]));
470  }
471  }
472  /* Issue a warning if extra factors of 1/z will be applied. */
473  if(numPoles<numZeros){
474  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
475  XLALPrintWarning("Filter has more zeros than poles\n");
476  XLALPrintWarning("\t%i poles added at complex origin\n",
477  numZeros-numPoles);
478  }
479  /* Issue a warning if any poles are outside |z|=1. */
480  for(i=0;i<numPoles;i++){
481  REAL8 zAbs=creal(poles[i])*creal(poles[i])+cimag(poles[i])*cimag(poles[i]);
482  if(zAbs>1.0){
483  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
484  XLALPrintWarning("Filter has pole outside of unit circle\n");
485  XLALPrintWarning("\tp_%i = %.8e + i*%.8e, |p_%i| = %.8e\n",i,
486  creal(poles[i]),cimag(poles[i]),i,zAbs);
487  }
488  }
489  }
490 #endif
491 
492  /* Everything seems okay, so initialize the filter. */
493  output=LALCalloc(1,sizeof(*output));
494  if ( ! output )
496  num = (numPoles>=numZeros) ? numZeros : numPoles;
497  numDirect+=num;
498  numRecurs+=num;
499 
500  output->deltaT=input->deltaT;
501  output->directCoef = XLALCreateREAL8Vector( numDirect );
502  output->recursCoef = XLALCreateREAL8Vector( numRecurs );
503  if ( ! output->directCoef || ! output->recursCoef )
504  {
505  XLALDestroyREAL8IIRFilter( output );
507  }
508  direct=output->directCoef->data;
509  recurs=output->recursCoef->data;
510 
511  /* Expand the denominator as a polynomial in z. */
512  *recurs=-1.0;
513  for(i=1;i<numRecurs;i++)
514  recurs[i]=0.0;
515  for(i=0;i<numPoles;i++){
516  INT4 j=numRecurs-1;
517  REAL8 x=creal(poles[i]);
518  REAL8 y=cimag(poles[i]);
519  if(y==0.0)
520  for(;j>0;j--)
521  recurs[j]-=x*recurs[j-1];
522  else if(y>0.0){
523  for(;j>1;j--)
524  recurs[j]-=2.0*x*recurs[j-1]-(x*x+y*y)*recurs[j-2];
525  recurs[j]-=2.0*x*recurs[j-1];
526  }
527  }
528 
529  /* Expand the numerator as a polynomial in z. */
530  for(i=0;i<numDirect;i++)
531  direct[i]=0.0;
532  direct[num-numZeros]=creal(input->gain);
533  for(i=0;i<numZeros;i++){
534  INT4 j=numDirect-1;
535  REAL8 x=creal(zeros[i]);
536  REAL8 y=cimag(zeros[i]);
537  if(y==0.0)
538  for(;j>num-numZeros;j--)
539  direct[j]-=x*direct[j-1];
540  else if(y>0.0){
541  for(;j>num-numZeros+1;j--)
542  direct[j]-=2.0*x*direct[j-1]-(x*x+y*y)*direct[j-2];
543  direct[j]-=2.0*x*direct[j-1];
544  }
545  }
546 
547  /* Initialize filter history. */
548  if(numDirect>=numRecurs)
549  num=numDirect-1;
550  else
551  num=numRecurs-1;
552  output->history = XLALCreateREAL8Vector( numRecurs );
553  if ( ! output->history )
554  {
555  XLALDestroyREAL8IIRFilter( output );
557  }
558  history=output->history->data;
559  for(i=0;i<num;i++)
560  history[i]=0.0;
561 
562  /* Normal exit */
563  return output;
564 }
565 
566 /** \see See \ref CreateIIRFilter_c for documentation */
568 {
569  COMPLEX16IIRFilter *output;
570  INT4 i; /* Index counter for zeros and poles. */
571  INT4 numZeros; /* The number of zeros. */
572  INT4 numPoles; /* The number of poles. */
573  INT4 numDirect; /* The number of direct filter coefficients. */
574  INT4 numRecurs; /* The number of recursive filter coefficients. */
575  INT4 num; /* An extra counter for error checking. */
576  COMPLEX16 *zeros; /* The zeros of the transfer function. */
577  COMPLEX16 *poles; /* The poles of the transfer function. */
578  REAL8 *direct; /* The direct filter coefficients. */
579  REAL8 *recurs; /* The recursive filter coefficients. */
580  COMPLEX16 *history; /* The filter history. */
581 
582  /* Make sure all the input structures have been initialized. */
583  if ( ! input )
585  if ( ! input->zeros || ! input->poles
586  || ! input->zeros->data || ! input->poles->data )
588 
589  numZeros=input->zeros->length;
590  numPoles=input->poles->length;
591  zeros=input->zeros->data;
592  poles=input->poles->data;
593 
594  /* Check that zeros are appropriately paired. Also, keep track of
595  the number of zeros at z=0, since these will reduce the number of
596  direct coefficients required. */
597  numDirect=1;
598  for(i=0,num=0;i<numZeros;i++)
599  if(cimag(zeros[i])==0.0){
600  num+=1;
601  if(creal(zeros[i])==0.0)
602  numDirect-=1;
603  }
604  else if(cimag(zeros[i])>0.0){
605  num+=2;
606 
607 #ifndef NDEBUG
609  /* Check to see that another zero is an actual conjugate.
610  This is not a foolproof test, as it can be fooled by
611  multiple zeros at the same location. */
612  INT4 j=0;
613  INT4 k=0;
614  REAL8 x = creal(zeros[i]) - creal(zeros[0]);
615  REAL8 y = cimag(zeros[i]) + cimag(zeros[0]);
616  REAL8 sep = x*x + y*y;
617  for(j=1;j<numZeros;j++){
618  x=creal(zeros[i])-creal(zeros[j]);
619  y=cimag(zeros[i])+cimag(zeros[j]);
620  if(sep>x*x+y*y){
621  sep=x*x+y*y;
622  k=j;
623  }
624  }
625  if(sep>LAL_REAL8_EPS){
626  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
627  XLALPrintWarning("Complex zero has no conjugate pair\n");
628  XLALPrintWarning("\tUnmatched zero z_%i = %.8e + i*%.8e\n",i,
629  creal(zeros[i]),cimag(zeros[i]));
630  XLALPrintWarning("\tNearest pair z_%i = %.8e + i*%.8e\n",k,
631  creal(zeros[k]),cimag(zeros[k]));
632  }
633  }
634 #endif
635  }
636 
637  if ( num != numZeros )
639 
640  /* Check that poles are appropriately paired. Also, keep track of
641  the number of poles at z=0, since these will reduce the number of
642  recursive coefficients required. */
643  numRecurs=1;
644  for(i=0,num=0;i<numPoles;i++)
645  if(cimag(poles[i])==0.0){
646  num+=1;
647  if(creal(poles[i])==0.0)
648  numRecurs-=1;
649  }
650  else if(cimag(poles[i])>0.0){
651  num+=2;
652 
653 #ifndef NDEBUG
655  /* Check to see that another pole is an actual conjugate.
656  This is not a foolproof test, as it can be fooled by
657  multiple poles at the same location. */
658  INT4 j=0;
659  INT4 k=0;
660  REAL8 x = creal(poles[i]) - creal(poles[0]);
661  REAL8 y = cimag(poles[i]) + cimag(poles[0]);
662  REAL8 sep = x*x + y*y;
663  for(j=1;j<numPoles;j++){
664  x=creal(poles[i])-creal(poles[j]);
665  y=cimag(poles[i])+cimag(poles[j]);
666  if(sep>x*x+y*y){
667  sep=x*x+y*y;
668  k=j;
669  }
670  }
671  if(sep>LAL_REAL8_EPS){
672  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
673  XLALPrintWarning("Complex pole has no conjugate pair\n");
674  XLALPrintWarning("\tUnmatched pole p_%i = %.8e + i*%.8e\n",i,
675  creal(poles[i]),cimag(poles[i]));
676  XLALPrintWarning("\tNearest pair p_%i = %.8e + i*%.8e\n",k,
677  creal(poles[k]),cimag(poles[k]));
678  }
679  }
680 #endif
681  }
682 
683  if ( num != numPoles )
685 
686 #ifndef NDEBUG
688  /* Issue a warning if the gain is nonreal. */
689  if(fabs(cimag(input->gain))>fabs(LAL_REAL8_EPS*creal(input->gain))){
690  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
691  XLALPrintWarning("Gain is non-real\n");
692  XLALPrintWarning("\tg = %.8e + i*%.8e\n", creal(input->gain), cimag(input->gain));
693  }
694  /* Issue a warning if there are any ``removeable'' poles. */
695  for(i=0;i<numPoles;i++){
696  INT4 j=0;
697  for(;j<numZeros;j++)
698  if((creal(poles[i])==creal(zeros[j]))&&(cimag(poles[i])==cimag(zeros[j]))){
699  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
700  XLALPrintWarning("Removeable pole\n");
701  XLALPrintWarning("\tp_%i = z_%i = %.8e + i*%.8e\n",i,j,
702  creal(poles[i]),cimag(poles[i]));
703  }
704  }
705  /* Issue a warning if extra factors of 1/z will be applied. */
706  if(numPoles<numZeros){
707  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
708  XLALPrintWarning("Filter has more zeros than poles\n");
709  XLALPrintWarning("\t%i poles added at complex origin\n",
710  numZeros-numPoles);
711  }
712  /* Issue a warning if any poles are outside |z|=1. */
713  for(i=0;i<numPoles;i++){
714  REAL8 zAbs=creal(poles[i])*creal(poles[i])+cimag(poles[i])*cimag(poles[i]);
715  if(zAbs>1.0){
716  XLALPrintWarning( "XLAL Warning - %s: ", __func__ );
717  XLALPrintWarning("Filter has pole outside of unit circle\n");
718  XLALPrintWarning("\tp_%i = %.8e + i*%.8e, |p_%i| = %.8e\n",i,
719  creal(poles[i]),cimag(poles[i]),i,zAbs);
720  }
721  }
722  }
723 #endif
724 
725  /* Everything seems okay, so initialize the filter. */
726  output=LALCalloc(1,sizeof(*output));
727  if ( ! output )
729  num = (numPoles>=numZeros) ? numZeros : numPoles;
730  numDirect+=num;
731  numRecurs+=num;
732 
733  output->deltaT=input->deltaT;
734  output->directCoef = XLALCreateREAL8Vector( numDirect );
735  output->recursCoef = XLALCreateREAL8Vector( numRecurs );
736  if ( ! output->directCoef || ! output->recursCoef )
737  {
740  }
741  direct=output->directCoef->data;
742  recurs=output->recursCoef->data;
743 
744  /* Expand the denominator as a polynomial in z. */
745  *recurs=-1.0;
746  for(i=1;i<numRecurs;i++)
747  recurs[i]=0.0;
748  for(i=0;i<numPoles;i++){
749  INT4 j=numRecurs-1;
750  REAL8 x=creal(poles[i]);
751  REAL8 y=cimag(poles[i]);
752  if(y==0.0)
753  for(;j>0;j--)
754  recurs[j]-=x*recurs[j-1];
755  else if(y>0.0){
756  for(;j>1;j--)
757  recurs[j]-=2.0*x*recurs[j-1]-(x*x+y*y)*recurs[j-2];
758  recurs[j]-=2.0*x*recurs[j-1];
759  }
760  }
761 
762  /* Expand the numerator as a polynomial in z. */
763  for(i=0;i<numDirect;i++)
764  direct[i]=0.0;
765  direct[num-numZeros]=creal(input->gain);
766  for(i=0;i<numZeros;i++){
767  INT4 j=numDirect-1;
768  REAL8 x=creal(zeros[i]);
769  REAL8 y=cimag(zeros[i]);
770  if(y==0.0)
771  for(;j>num-numZeros;j--)
772  direct[j]-=x*direct[j-1];
773  else if(y>0.0){
774  for(;j>num-numZeros+1;j--)
775  direct[j]-=2.0*x*direct[j-1]-(x*x+y*y)*direct[j-2];
776  direct[j]-=2.0*x*direct[j-1];
777  }
778  }
779 
780  /* Initialize filter history. */
781  if(numDirect>=numRecurs)
782  num=numDirect-1;
783  else
784  num=numRecurs-1;
785  output->history = XLALCreateCOMPLEX16Vector( numRecurs );
786  if ( ! output->history )
787  {
790  }
791  history=output->history->data;
792  for(i=0;i<num;i++)
793  history[i]=0.0;
794 
795  /* Normal exit */
796  return output;
797 }
798 
799 
800 /**
801  * Deprecated.
802  * \deprecated Use XLALCreateREAL4IIRFilter() instead
803  */
805  REAL4IIRFilter **output,
806  COMPLEX8ZPGFilter *input )
807 {
808  INITSTATUS(stat);
809 
810  /* Make sure all the input structures have been initialized. */
811  ASSERT(input,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
812  ASSERT(input->zeros,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
813  ASSERT(input->zeros->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
814  ASSERT(input->poles,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
815  ASSERT(input->poles->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
816 
817  /* Make sure that the output handle exists, but points to a null
818  pointer. */
819  ASSERT(output,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
820  ASSERT(!*output,stat,IIRFILTERH_EOUT,IIRFILTERH_MSGEOUT);
821 
822  *output = XLALCreateREAL4IIRFilter( input );
823  if ( ! *output )
824  {
825  int code = xlalErrno & ~XLAL_EFUNC;
826  XLALClearErrno();
827  switch ( code )
828  {
829  case XLAL_EINVAL:
830  ABORT(stat,IIRFILTERH_EPAIR,IIRFILTERH_MSGEPAIR);
831  case XLAL_ENOMEM:
832  ABORT(stat,IIRFILTERH_EMEM,IIRFILTERH_MSGEMEM);
833  default:
834  ABORTXLAL(stat);
835  }
836  }
837 
838  RETURN(stat);
839 }
840 
841 
842 /**
843  * Deprecated.
844  * \deprecated Use XLALCreateREAL8IIRFilter() instead
845  */
847  REAL8IIRFilter **output,
848  COMPLEX16ZPGFilter *input )
849 {
850  INITSTATUS(stat);
851 
852  /* Make sure all the input structures have been initialized. */
853  ASSERT(input,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
854  ASSERT(input->zeros,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
855  ASSERT(input->zeros->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
856  ASSERT(input->poles,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
857  ASSERT(input->poles->data,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
858 
859  /* Make sure that the output handle exists, but points to a null
860  pointer. */
861  ASSERT(output,stat,IIRFILTERH_ENUL,IIRFILTERH_MSGENUL);
862  ASSERT(!*output,stat,IIRFILTERH_EOUT,IIRFILTERH_MSGEOUT);
863 
864  *output = XLALCreateREAL8IIRFilter( input );
865  if ( ! *output )
866  {
867  int code = xlalErrno & ~XLAL_EFUNC;
868  XLALClearErrno();
869  switch ( code )
870  {
871  case XLAL_EINVAL:
872  ABORT(stat,IIRFILTERH_EPAIR,IIRFILTERH_MSGEPAIR);
873  case XLAL_ENOMEM:
874  ABORT(stat,IIRFILTERH_EMEM,IIRFILTERH_MSGEMEM);
875  default:
876  ABORTXLAL(stat);
877  }
878  }
879 
880  RETURN(stat);
881 }
882 /** @} */
This structure stores the direct and recursive REAL8 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:168
REAL4IIRFilter * XLALCreateREAL4IIRFilter(COMPLEX8ZPGFilter *input)
#define LAL_REAL4_EPS
Difference between 1 and the next resolvable REAL4 2^-23.
Definition: LALConstants.h:57
REAL8Vector * directCoef
The direct filter coefficients.
Definition: IIRFilter.h:171
int XLALPrintWarning(const char *fmt,...)
Definition: XLALError.c:79
REAL8Vector * recursCoef
The recursive filter coefficients.
Definition: IIRFilter.h:172
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
REAL8Vector * directCoef
The direct filter coefficients.
Definition: IIRFilter.h:187
#define IIRFILTERH_EPAIR
Input has unpaired nonreal poles or zeros.
Definition: IIRFilter.h:134
#define xlalErrno
Modifiable lvalue containing the XLAL error number.
Definition: XLALError.h:572
REAL8 * data
Pointer to the data array.
Definition: LALDatatypes.h:159
COMPLEX8Vector * zeros
Definition: LALDatatypes.h:924
void XLALDestroyREAL4IIRFilter(REAL4IIRFilter *filter)
#define IIRFILTERH_EOUT
Output handle points to a non-null pointer.
Definition: IIRFilter.h:132
COMPLEX16Vector * history
The previous values of w.
Definition: IIRFilter.h:189
This structure stores the direct and recursive REAL4 filter coefficients, as well as the history of t...
Definition: IIRFilter.h:152
void XLALDestroyREAL8IIRFilter(REAL8IIRFilter *filter)
void XLALDestroyCOMPLEX16IIRFilter(COMPLEX16IIRFilter *filter)
REAL4Vector * directCoef
The direct filter coefficients.
Definition: IIRFilter.h:155
float complex COMPLEX8
Single-precision floating-point complex number (8 bytes total)
See DATATYPE-ZPGFilter types for details.
Definition: LALDatatypes.h:930
REAL8Vector * history
The previous values of w.
Definition: IIRFilter.h:173
UINT4 length
Number of elements in array.
Definition: LALDatatypes.h:176
See DATATYPE-ZPGFilter types for details.
Definition: LALDatatypes.h:921
REAL4Vector * history
The previous values of w.
Definition: IIRFilter.h:157
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
float REAL4
Single precision real floating-point number (4 bytes).
#define IIRFILTERH_ENUL
Unexpected null pointer in arguments.
Definition: IIRFilter.h:131
This structure stores the direct and recursive REAL8 filter coefficients, as well as the complex-valu...
Definition: IIRFilter.h:184
REAL8Vector * recursCoef
The recursive filter coefficients.
Definition: IIRFilter.h:188
UINT4 length
Number of elements in array.
Definition: LALDatatypes.h:167
#define LAL_REAL8_EPS
Difference between 1 and the next resolvable REAL8 2^-52.
Definition: LALConstants.h:61
#define RETURN(statusptr)
#define INITSTATUS(statusptr)
int XLALClearErrno(void)
Clears the XLAL error number, returns the old value.
Definition: XLALError.c:363
Invalid pointer.
Definition: XLALError.h:409
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
REAL4Vector * XLALCreateREAL4Vector(UINT4 length)
#define XLAL_ERROR_NULL(...)
Macro to invoke a failure from a XLAL routine returning a pointer.
Definition: XLALError.h:714
COMPLEX16Vector * zeros
Definition: LALDatatypes.h:933
void LALCreateREAL8IIRFilter(LALStatus *stat, REAL8IIRFilter **output, COMPLEX16ZPGFilter *input)
Deprecated.
double REAL8
Double precision real floating-point number (8 bytes).
#define ABORT(statusptr, code, mesg)
COMPLEX16 * data
Pointer to the data array.
Definition: LALDatatypes.h:177
#define lalDebugLevel
Definition: LALDebugLevel.h:56
REAL4 * data
Pointer to the data array.
Definition: LALDatatypes.h:150
COMPLEX16IIRFilter * XLALCreateCOMPLEX16IIRFilter(COMPLEX16ZPGFilter *input)
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
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
#define LALCalloc(m, n)
Definition: LALMalloc.h:94
COMPLEX16Vector * poles
Definition: LALDatatypes.h:934
#define IIRFILTERH_EMEM
Memory allocation error.
Definition: IIRFilter.h:133
REAL8IIRFilter * XLALCreateREAL8IIRFilter(COMPLEX16ZPGFilter *input)
double complex COMPLEX16
Double-precision floating-point complex number (16 bytes total)
COMPLEX8Vector * poles
Definition: LALDatatypes.h:925
#define ABORTXLAL(sp)
Memory allocation error.
Definition: XLALError.h:408
COMPLEX16Vector * XLALCreateCOMPLEX16Vector(UINT4 length)
int32_t INT4
Four-byte signed integer.
Internal function call failed bit: "or" this with existing error number.
Definition: XLALError.h:463
Invalid argument.
Definition: XLALError.h:410
void LALCreateREAL4IIRFilter(LALStatus *stat, REAL4IIRFilter **output, COMPLEX8ZPGFilter *input)
Deprecated.
#define ASSERT(assertion, statusptr, code, mesg)
COMPLEX8 * data
Pointer to the data array.
Definition: LALDatatypes.h:168
enable warning messages
Definition: LALDebugLevel.h:45
REAL4Vector * recursCoef
The recursive filter coefficients.
Definition: IIRFilter.h:156