openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
iir.c
Go to the documentation of this file.
1 /*
2  * COPYRIGHT
3  *
4  * liir - Recursive digital filter functions
5  * Copyright (C) 2007 Exstrom Laboratories LLC
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * A copy of the GNU General Public License is available on the internet at:
18  *
19  * http://www.gnu.org/copyleft/gpl.html
20  *
21  * or you can write to:
22  *
23  * The Free Software Foundation, Inc.
24  * 675 Mass Ave
25  * Cambridge, MA 02139, USA
26  *
27  * You can contact Exstrom Laboratories LLC via Email at:
28  *
29  * stefan(AT)exstrom.com
30  *
31  * or you can write to:
32  *
33  * Exstrom Laboratories LLC
34  * P.O. Box 7651
35  * Longmont, CO 80501, USA
36  *
37  */
38 
39 #include <stdlib.h>
40 #include <stdio.h>
41 #include <string.h>
42 #ifndef __USE_GNU
43 #define __USE_GNU
44 #endif
45 #include <math.h>
46 #include "iir.h"
47 
48 #ifndef M_PIl
49 #define M_PIl (3.14159265358979323846264338327950288)
50 #endif
51 
52 
53 /**********************************************************************
54  binomial_mult - multiplies a series of binomials together and returns
55  the coefficients of the resulting polynomial.
56 
57  The multiplication has the following form:
58 
59  (x+p[0])*(x+p[1])*...*(x+p[n-1])
60 
61  The p[i] coefficients are assumed to be complex and are passed to the
62  function as a pointer to an array of long double s of length 2n.
63 
64  The resulting polynomial has the following form:
65 
66  x^n + a[0]*x^n-1 + a[1]*x^n-2 + ... +a[n-2]*x + a[n-1]
67 
68  The a[i] coefficients can in general be complex but should in most
69  cases turn out to be real. The a[i] coefficients are returned by the
70  function as a pointer to an array of long double s of length 2n. Storage
71  for the array is allocated by the function and should be freed by the
72  calling program when no longer needed.
73 
74  Function arguments:
75 
76  n - The number of binomials to multiply
77  p - Pointer to an array of long double s where p[2i] (i=0...n-1) is
78  assumed to be the real part of the coefficient of the ith binomial
79  and p[2i+1] is assumed to be the imaginary part. The overall size
80  of the array is then 2n.
81 */
82 
83 long double *binomial_mult( int n, long double *p )
84 {
85  int i, j;
86  long double *a;
87 
88  a = (long double *)calloc( 2 * n, sizeof(long double ) );
89  if( a == NULL ) return( NULL );
90 
91  for( i = 0; i < n; ++i )
92  {
93  for( j = i; j > 0; --j )
94  {
95  a[2*j] += p[2*i] * a[2*(j-1)] - p[2*i+1] * a[2*(j-1)+1];
96  a[2*j+1] += p[2*i] * a[2*(j-1)+1] + p[2*i+1] * a[2*(j-1)];
97  }
98  a[0] += p[2*i];
99  a[1] += p[2*i+1];
100  }
101  return( a );
102 }
103 
104 
105 /**********************************************************************
106  trinomial_mult - multiplies a series of trinomials together and returns
107  the coefficients of the resulting polynomial.
108 
109  The multiplication has the following form:
110 
111  (x^2 + b[0]x + c[0])*(x^2 + b[1]x + c[1])*...*(x^2 + b[n-1]x + c[n-1])
112 
113  The b[i] and c[i] coefficients are assumed to be complex and are passed
114  to the function as a pointers to arrays of long double s of length 2n. The real
115  part of the coefficients are stored in the even numbered elements of the
116  array and the imaginary parts are stored in the odd numbered elements.
117 
118  The resulting polynomial has the following form:
119 
120  x^2n + a[0]*x^2n-1 + a[1]*x^2n-2 + ... +a[2n-2]*x + a[2n-1]
121 
122  The a[i] coefficients can in general be complex but should in most cases
123  turn out to be real. The a[i] coefficients are returned by the function as
124  a pointer to an array of long double s of length 4n. The real and imaginary
125  parts are stored, respectively, in the even and odd elements of the array.
126  Storage for the array is allocated by the function and should be freed by
127  the calling program when no longer needed.
128 
129  Function arguments:
130 
131  n - The number of trinomials to multiply
132  b - Pointer to an array of long double s of length 2n.
133  c - Pointer to an array of long double s of length 2n.
134 */
135 
136 long double *trinomial_mult( int n, long double *b, long double *c )
137 {
138  int i, j;
139  long double *a;
140 
141  a = (long double *)calloc( 4 * n, sizeof(long double ) );
142  if( a == NULL ) return( NULL );
143 
144  a[2] = c[0];
145  a[3] = c[1];
146  a[0] = b[0];
147  a[1] = b[1];
148 
149  for( i = 1; i < n; ++i )
150  {
151  a[2*(2*i+1)] += c[2*i]*a[2*(2*i-1)] - c[2*i+1]*a[2*(2*i-1)+1];
152  a[2*(2*i+1)+1] += c[2*i]*a[2*(2*i-1)+1] + c[2*i+1]*a[2*(2*i-1)];
153 
154  for( j = 2*i; j > 1; --j )
155  {
156  a[2*j] += b[2*i] * a[2*(j-1)] - b[2*i+1] * a[2*(j-1)+1] +
157  c[2*i] * a[2*(j-2)] - c[2*i+1] * a[2*(j-2)+1];
158  a[2*j+1] += b[2*i] * a[2*(j-1)+1] + b[2*i+1] * a[2*(j-1)] +
159  c[2*i] * a[2*(j-2)+1] + c[2*i+1] * a[2*(j-2)];
160  }
161 
162  a[2] += b[2*i] * a[0] - b[2*i+1] * a[1] + c[2*i];
163  a[3] += b[2*i] * a[1] + b[2*i+1] * a[0] + c[2*i+1];
164  a[0] += b[2*i];
165  a[1] += b[2*i+1];
166  }
167 
168  return( a );
169 }
170 
171 
172 /**********************************************************************
173  dcof_bwlp - calculates the d coefficients for a butterworth lowpass
174  filter. The coefficients are returned as an array of long double s.
175 
176 */
177 
178 long double *dcof_bwlp( int n, long double fcf )
179 {
180  int k; // loop variables
181  long double theta; // M_PIl * fcf / 2.0
182  long double st; // sine of theta
183  long double ct; // cosine of theta
184  long double parg; // pole angle
185  long double sparg; // sine of the pole angle
186  long double cparg; // cosine of the pole angle
187  long double a; // workspace variable
188  long double *rcof; // binomial coefficients
189  long double *dcof; // dk coefficients
190 
191  rcof = (long double *)calloc( 2 * n, sizeof(long double ) );
192  if( rcof == NULL ) return( NULL );
193 
194  theta = M_PIl * fcf;
195  st = sinl(theta);
196  ct = cosl(theta);
197 
198  for( k = 0; k < n; ++k )
199  {
200  parg = M_PIl * (long double )(2*k+1)/(long double )(2*n);
201  sparg = sinl(parg);
202  cparg = cosl(parg);
203  a = 1.0 + st*sparg;
204  rcof[2*k] = -ct/a;
205  rcof[2*k+1] = -st*cparg/a;
206  }
207 
208  dcof = binomial_mult( n, rcof );
209  free( rcof );
210 
211  dcof[1] = dcof[0];
212  dcof[0] = 1.0;
213  for( k = 3; k <= n; ++k )
214  dcof[k] = dcof[2*k-2];
215  return( dcof );
216 }
217 
218 /**********************************************************************
219  dcof_bwhp - calculates the d coefficients for a butterworth highpass
220  filter. The coefficients are returned as an array of long double s.
221 
222 */
223 
224 long double *dcof_bwhp( int n, long double fcf )
225 {
226  return( dcof_bwlp( n, fcf ) );
227 }
228 
229 
230 /**********************************************************************
231  dcof_bwbp - calculates the d coefficients for a butterworth bandpass
232  filter. The coefficients are returned as an array of long double s.
233 
234 */
235 
236 long double *dcof_bwbp( int n, long double f1f, long double f2f )
237 {
238  int k; // loop variables
239  long double theta; // M_PIl * (f2f - f1f) / 2.0
240  long double cp; // cosine of phi
241  long double st; // sine of theta
242  long double ct; // cosine of theta
243  long double s2t; // sine of 2*theta
244  long double c2t; // cosine 0f 2*theta
245  long double *rcof; // z^-2 coefficients
246  long double *tcof; // z^-1 coefficients
247  long double *dcof; // dk coefficients
248  long double parg; // pole angle
249  long double sparg; // sine of pole angle
250  long double cparg; // cosine of pole angle
251  long double a; // workspace variables
252 
253  cp = cosl(M_PIl * (f2f + f1f) / 2.0);
254  theta = M_PIl * (f2f - f1f) / 2.0;
255  st = sinl(theta);
256  ct = cosl(theta);
257  s2t = 2.0*st*ct; // sine of 2*theta
258  c2t = 2.0*ct*ct - 1.0; // cosine of 2*theta
259 
260  rcof = (long double *)calloc( 2 * n, sizeof(long double ) );
261  tcof = (long double *)calloc( 2 * n, sizeof(long double ) );
262 
263  for( k = 0; k < n; ++k )
264  {
265  parg = M_PIl * (long double )(2*k+1)/(long double )(2*n);
266  sparg = sinl(parg);
267  cparg = cosl(parg);
268  a = 1.0 + s2t*sparg;
269  rcof[2*k] = c2t/a;
270  rcof[2*k+1] = s2t*cparg/a;
271  tcof[2*k] = -2.0*cp*(ct+st*sparg)/a;
272  tcof[2*k+1] = -2.0*cp*st*cparg/a;
273  }
274 
275  dcof = trinomial_mult( n, tcof, rcof );
276  free( tcof );
277  free( rcof );
278 
279  dcof[1] = dcof[0];
280  dcof[0] = 1.0;
281  for( k = 3; k <= 2*n; ++k )
282  dcof[k] = dcof[2*k-2];
283  return( dcof );
284 }
285 
286 /**********************************************************************
287  dcof_bwbs - calculates the d coefficients for a butterworth bandstop
288  filter. The coefficients are returned as an array of long double s.
289 
290 */
291 
292 long double *dcof_bwbs( int n, long double f1f, long double f2f )
293 {
294  int k; // loop variables
295  long double theta; // M_PIl * (f2f - f1f) / 2.0
296  long double cp; // cosine of phi
297  long double st; // sine of theta
298  long double ct; // cosine of theta
299  long double s2t; // sine of 2*theta
300  long double c2t; // cosine 0f 2*theta
301  long double *rcof; // z^-2 coefficients
302  long double *tcof; // z^-1 coefficients
303  long double *dcof; // dk coefficients
304  long double parg; // pole angle
305  long double sparg; // sine of pole angle
306  long double cparg; // cosine of pole angle
307  long double a; // workspace variables
308 
309  cp = cosl(M_PIl * (f2f + f1f) / 2.0);
310  theta = M_PIl * (f2f - f1f) / 2.0;
311  st = sinl(theta);
312  ct = cosl(theta);
313  s2t = 2.0*st*ct; // sine of 2*theta
314  c2t = 2.0*ct*ct - 1.0; // cosine 0f 2*theta
315 
316  rcof = (long double *)calloc( 2 * n, sizeof(long double ) );
317  tcof = (long double *)calloc( 2 * n, sizeof(long double ) );
318 
319  for( k = 0; k < n; ++k )
320  {
321  parg = M_PIl * (long double )(2*k+1)/(long double )(2*n);
322  sparg = sinl(parg);
323  cparg = cosl(parg);
324  a = 1.0 + s2t*sparg;
325  rcof[2*k] = c2t/a;
326  rcof[2*k+1] = -s2t*cparg/a;
327  tcof[2*k] = -2.0*cp*(ct+st*sparg)/a;
328  tcof[2*k+1] = 2.0*cp*st*cparg/a;
329  }
330 
331  dcof = trinomial_mult( n, tcof, rcof );
332  free( tcof );
333  free( rcof );
334 
335  dcof[1] = dcof[0];
336  dcof[0] = 1.0;
337  for( k = 3; k <= 2*n; ++k )
338  dcof[k] = dcof[2*k-2];
339  return( dcof );
340 }
341 
342 /**********************************************************************
343  ccof_bwlp - calculates the c coefficients for a butterworth lowpass
344  filter. The coefficients are returned as an array of integers.
345 
346 */
347 
348 int *ccof_bwlp( int n )
349 {
350  int *ccof;
351  int m;
352  int i;
353 
354  ccof = (int *)calloc( n+1, sizeof(int) );
355  if( ccof == NULL ) return( NULL );
356 
357  ccof[0] = 1;
358  ccof[1] = n;
359  m = n/2;
360  for( i=2; i <= m; ++i)
361  {
362  ccof[i] = (n-i+1)*ccof[i-1]/i;
363  ccof[n-i]= ccof[i];
364  }
365  ccof[n-1] = n;
366  ccof[n] = 1;
367 
368  return( ccof );
369 }
370 
371 /**********************************************************************
372  ccof_bwhp - calculates the c coefficients for a butterworth highpass
373  filter. The coefficients are returned as an array of integers.
374 
375 */
376 
377 int *ccof_bwhp( int n )
378 {
379  int *ccof;
380  int i;
381 
382  ccof = ccof_bwlp( n );
383  if( ccof == NULL ) return( NULL );
384 
385  for( i = 0; i <= n; ++i)
386  if( i % 2 ) ccof[i] = -ccof[i];
387 
388  return( ccof );
389 }
390 
391 /**********************************************************************
392  ccof_bwbp - calculates the c coefficients for a butterworth bandpass
393  filter. The coefficients are returned as an array of integers.
394 
395 */
396 
397 int *ccof_bwbp( int n )
398 {
399  int *tcof;
400  int *ccof;
401  int i;
402 
403  ccof = (int *)calloc( 2*n+1, sizeof(int) );
404  if( ccof == NULL ) return( NULL );
405 
406  tcof = ccof_bwhp(n);
407  if( tcof == NULL ) return( NULL );
408 
409  for( i = 0; i < n; ++i)
410  {
411  ccof[2*i] = tcof[i];
412  ccof[2*i+1] = 0.0;
413  }
414  ccof[2*n] = tcof[n];
415 
416  free( tcof );
417  return( ccof );
418 }
419 
420 /**********************************************************************
421  ccof_bwbs - calculates the c coefficients for a butterworth bandstop
422  filter. The coefficients are returned as an array of integers.
423 
424 */
425 
426 long double *ccof_bwbs( int n, long double f1f, long double f2f )
427 {
428  long double alpha;
429  long double *ccof;
430  int i, j;
431 
432  alpha = -2.0 * cosl(M_PIl * (f2f + f1f) / 2.0) / cos(M_PIl * (f2f - f1f) / 2.0);
433 
434  ccof = (long double *)calloc( 2*n+1, sizeof(long double ) );
435 
436  ccof[0] = 1.0;
437 
438  ccof[2] = 1.0;
439  ccof[1] = alpha;
440 
441  for( i = 1; i < n; ++i )
442  {
443  ccof[2*i+2] += ccof[2*i];
444  for( j = 2*i; j > 1; --j )
445  ccof[j+1] += alpha * ccof[j] + ccof[j-1];
446 
447  ccof[2] += alpha * ccof[1] + 1.0;
448  ccof[1] += alpha;
449  }
450 
451  return( ccof );
452 }
453 
454 /**********************************************************************
455  sf_bwlp - calculates the scaling factor for a butterworth lowpass filter.
456  The scaling factor is what the c coefficients must be multiplied by so
457  that the filter response has a maximum value of 1.
458 
459 */
460 
461 long double sf_bwlp( int n, long double fcf )
462 {
463  int m, k; // loop variables
464  long double omega; // M_PIl * fcf
465  long double fomega; // function of omega
466  long double parg0; // zeroth pole angle
467  long double sf; // scaling factor
468 
469  omega = M_PIl * fcf;
470  fomega = sinl(omega);
471  parg0 = M_PIl / (long double )(2*n);
472 
473  m = n / 2;
474  sf = 1.0;
475  for( k = 0; k < n/2; ++k )
476  sf *= 1.0 + fomega * sinl((long double )(2*k+1)*parg0);
477 
478  fomega = sinl(omega / 2.0);
479 
480  if( n % 2 ) sf *= fomega + cosl(omega / 2.0);
481  sf = pow( fomega, n ) / sf;
482 
483  return(sf);
484 }
485 
486 /**********************************************************************
487  sf_bwhp - calculates the scaling factor for a butterworth highpass filter.
488  The scaling factor is what the c coefficients must be multiplied by so
489  that the filter response has a maximum value of 1.
490 
491 */
492 
493 long double sf_bwhp( int n, long double fcf )
494 {
495  int m, k; // loop variables
496  long double omega; // M_PIl * fcf
497  long double fomega; // function of omega
498  long double parg0; // zeroth pole angle
499  long double sf; // scaling factor
500 
501  omega = M_PIl * fcf;
502  fomega = sinl(omega);
503  parg0 = M_PIl / (long double )(2*n);
504 
505  m = n / 2;
506  sf = 1.0;
507  for( k = 0; k < n/2; ++k )
508  sf *= 1.0 + fomega * sinl((long double )(2*k+1)*parg0);
509 
510  fomega = cosl(omega / 2.0);
511 
512  if( n % 2 ) sf *= fomega + sinl(omega / 2.0);
513  sf = pow( fomega, n ) / sf;
514 
515  return(sf);
516 }
517 
518 /**********************************************************************
519  sf_bwbp - calculates the scaling factor for a butterworth bandpass filter.
520  The scaling factor is what the c coefficients must be multiplied by so
521  that the filter response has a maximum value of 1.
522 
523 */
524 
525 long double sf_bwbp( int n, long double f1f, long double f2f )
526 {
527  int k; // loop variables
528  long double ctt; // cotangent of theta
529  long double sfr, sfi; // real and imaginary parts of the scaling factor
530  long double parg; // pole angle
531  long double sparg; // sine of pole angle
532  long double cparg; // cosine of pole angle
533  long double a, b, c; // workspace variables
534 
535  ctt = 1.0 / tanl(M_PIl * (f2f - f1f) / 2.0);
536  sfr = 1.0;
537  sfi = 0.0;
538 
539  for( k = 0; k < n; ++k )
540  {
541  parg = M_PIl * (long double )(2*k+1)/(long double )(2*n);
542  sparg = ctt + sinl(parg);
543  cparg = cosl(parg);
544  a = (sfr + sfi)*(sparg - cparg);
545  b = sfr * sparg;
546  c = -sfi * cparg;
547  sfr = b - c;
548  sfi = a - b - c;
549  }
550 
551  return( 1.0 / sfr );
552 }
553 
554 /**********************************************************************
555  sf_bwbs - calculates the scaling factor for a butterworth bandstop filter.
556  The scaling factor is what the c coefficients must be multiplied by so
557  that the filter response has a maximum value of 1.
558 
559 */
560 
561 long double sf_bwbs( int n, long double f1f, long double f2f )
562 {
563  int k; // loop variables
564  long double tt; // tangent of theta
565  long double sfr, sfi; // real and imaginary parts of the scaling factor
566  long double parg; // pole angle
567  long double sparg; // sine of pole angle
568  long double cparg; // cosine of pole angle
569  long double a, b, c; // workspace variables
570 
571  tt = tanl(M_PIl * (f2f - f1f) / 2.0);
572  sfr = 1.0;
573  sfi = 0.0;
574 
575  for( k = 0; k < n; ++k )
576  {
577  parg = M_PIl * (long double )(2*k+1)/(long double )(2*n);
578  sparg = tt + sinl(parg);
579  cparg = cosl(parg);
580  a = (sfr + sfi)*(sparg - cparg);
581  b = sfr * sparg;
582  c = -sfi * cparg;
583  sfr = b - c;
584  sfi = a - b - c;
585  }
586 
587  return( 1.0 / sfr );
588 }
long double * trinomial_mult(int n, long double *b, long double *c)
Definition: iir.c:136
long double sf_bwhp(int n, long double fcf)
Definition: iir.c:493
long double * dcof_bwbp(int n, long double f1f, long double f2f)
Definition: iir.c:236
int * ccof_bwbp(int n)
Definition: iir.c:397
long double * binomial_mult(int n, long double *p)
Definition: iir.c:83
long double * dcof_bwbs(int n, long double f1f, long double f2f)
Definition: iir.c:292
int * ccof_bwhp(int n)
Definition: iir.c:377
long double sf_bwbs(int n, long double f1f, long double f2f)
Definition: iir.c:561
long double sf_bwbp(int n, long double f1f, long double f2f)
Definition: iir.c:525
#define M_PIl
Definition: iir.c:49
int * ccof_bwlp(int n)
Definition: iir.c:348
long double * dcof_bwlp(int n, long double fcf)
Definition: iir.c:178
long double * ccof_bwbs(int n, long double f1f, long double f2f)
Definition: iir.c:426
long double sf_bwlp(int n, long double fcf)
Definition: iir.c:461
long double * dcof_bwhp(int n, long double fcf)
Definition: iir.c:224