libdspl-2.0
Digital Signal Processing Algorithm Library
filter_iir.c
1 /*
2 * Copyright (c) 2015-2019 Sergey Bakhurin
3 * Digital Signal Processing Library [http://dsplib.org]
4 *
5 * This file is part of libdspl-2.0.
6 *
7 * is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * DSPL 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 * You should have received a copy of the GNU Lesser General Public License
18 * along with Foobar. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #include <stdlib.h>
22 #include <stdio.h>
23 #include <string.h>
24 #include "dspl.h"
25 #include "dspl_internal.h"
26 
27 
28 
29 #ifdef DOXYGEN_ENGLISH
30 
117 #endif
118 #ifdef DOXYGEN_RUSSIAN
119 
209 #endif
210 int DSPL_API bilinear(double* bs, double* as, int ord, double* bz, double* az)
211 {
212  double c[2] = {1.0, -1.0};
213  double d[2] = {1.0, 1.0};
214  return ratcompos(bs, as, ord, c, d, 1, bz, az);
215 }
216 
217 
218 
219 
220 
221 #ifdef DOXYGEN_ENGLISH
222 
310 #endif
311 #ifdef DOXYGEN_RUSSIAN
312 
403 #endif
404 int DSPL_API iir(double rp, double rs, int ord, double w0, double w1,
405  int type, double* b, double* a)
406 {
407  double *bs = NULL;
408  double *as = NULL;
409  double *bt = NULL;
410  double *at = NULL;
411  double wa0, wa1, ws;
412  int err, ord_ap = ord;
413  int i;
414 
415  if(((type & DSPL_FILTER_TYPE_MASK) == DSPL_FILTER_LPF) ||
416  ((type & DSPL_FILTER_TYPE_MASK) == DSPL_FILTER_HPF))
417  {
418  bs = (double*)malloc((ord_ap+1)*sizeof(double));
419  as = (double*)malloc((ord_ap+1)*sizeof(double));
420  bt = (double*)malloc((ord_ap+1)*sizeof(double));
421  at = (double*)malloc((ord_ap+1)*sizeof(double));
422  }
423 
424 
425  if(((type & DSPL_FILTER_TYPE_MASK) == DSPL_FILTER_BPASS) ||
426  ((type & DSPL_FILTER_TYPE_MASK) == DSPL_FILTER_BSTOP))
427  {
428  if(ord % 2)
429  return ERROR_FILTER_ORD_BP;
430  else
431  {
432  ord_ap = ord / 2;
433  bs = (double*)malloc((ord_ap + 1)*sizeof(double));
434  as = (double*)malloc((ord_ap + 1)*sizeof(double));
435  bt = (double*)malloc((ord + 1)*sizeof(double));
436  at = (double*)malloc((ord + 1)*sizeof(double));
437  }
438  }
439  err = iir_ap(rp, rs, ord_ap, type, bs, as);
440  if(err != RES_OK)
441  goto error_proc;
442 
443  /* frequency transformation */
444  wa0 = tan(w0 * M_PI * 0.5);
445  wa1 = tan(w1 * M_PI * 0.5);
446 
447  switch(type & DSPL_FILTER_TYPE_MASK)
448  {
449 
450  case DSPL_FILTER_LPF:
451  err = low2low(bs, as, ord_ap, 1.0, wa0, bt, at);
452  break;
453 
454  case DSPL_FILTER_HPF:
455  ws = filter_ws1(ord_ap, rp, rs, type);
456  err = low2low( bs, as, ord_ap, 1.0, 1.0 / ws, bs, as);
457  err = low2high(bs, as, ord_ap, 1.0, wa0, bt, at);
458  break;
459 
460  case DSPL_FILTER_BPASS:
461  err = low2bp(bs, as, ord_ap, 1.0, wa0, wa1, bt, at);
462  break;
463 
464  case DSPL_FILTER_BSTOP:
465  /* need frequency transform ws -> 1 rad/s */
466 
467  ws = filter_ws1(ord_ap, rp, rs, type);
468  err = low2low( bs, as, ord_ap, 1.0, 1.0 / ws, bs, as);
469  err = low2bs(bs, as, ord_ap, 1.0, wa0, wa1, bt, at);
470  break;
471 
472  default:
473  err = ERROR_FILTER_TYPE;
474  break;
475  }
476  if(err != RES_OK)
477  goto error_proc;
478 
479 
480  err = bilinear(bt, at, ord, b, a);
481 
482  for(i = 1; i <= ord; i++)
483  {
484  a[i] /= a[0];
485  b[i] /= a[0];
486  }
487  b[0] /= a[0];
488  a[0] = 1.0;
489 
490 error_proc:
491 
492  if(bs)
493  free(bs);
494  if(as)
495  free(as);
496  if(bt)
497  free(bt);
498  if(at)
499  free(at);
500 
501  return err;
502 
503 }
504 
505 
506 
507 
508 
509 
510 /******************************************************************************
511 Analog prototype for IIR
512 *******************************************************************************/
513 int iir_ap(double rp, double rs, int ord, int type, double* b, double* a)
514 {
515  int err;
516  switch(type & DSPL_FILTER_APPROX_MASK)
517  {
518  case DSPL_FILTER_BUTTER:
519  err = butter_ap(rp, ord, b, a);
520  break;
521  case DSPL_FILTER_CHEBY1:
522  err = cheby1_ap(rp, ord, b, a);
523  break;
524  case DSPL_FILTER_CHEBY2:
525  err = cheby2_ap_wp1(rp, rs, ord, b, a);
526  break;
527  case DSPL_FILTER_ELLIP:
528  err = ellip_ap(rp, rs, ord, b, a);
529  break;
530  default:
531  err = ERROR_FILTER_APPROX;
532  }
533  return err;
534 }
535 
536 
537 
538 
int DSPL_API cheby1_ap(double rp, int ord, double *b, double *a)
Function calculates the transfer function coefficients of analog normalized lowpass Chebyshev type 1...
Definition: filter_ap.c:601
int DSPL_API ratcompos(double *b, double *a, int n, double *c, double *d, int p, double *beta, double *alpha)
Rational composition.
Definition: filter_ft.c:602
int DSPL_API low2low(double *b, double *a, int ord, double w0, double w1, double *beta, double *alpha)
Definition: filter_ft.c:426
int DSPL_API bilinear(double *bs, double *as, int ord, double *bz, double *az)
Transform a s-plane analog filter transfer function to the digital filter transfer function .
Definition: filter_iir.c:210
int DSPL_API iir(double rp, double rs, int ord, double w0, double w1, int type, double *b, double *a)
Digital IIR filter design.
Definition: filter_iir.c:404
#define ERROR_FILTER_TYPE
Unknown filter type. The library supports the following filter types: lowpass, highpass,...
Definition: dspl.h:518
int DSPL_API low2high(double *b, double *a, int ord, double w0, double w1, double *beta, double *alpha)
Lowpass to highpass filter frequency transform.
Definition: filter_ft.c:280
#define ERROR_FILTER_APPROX
Unknown approximation type of digital or analog filter. This error occurs when masks such as a digita...
Definition: dspl.h:512
int DSPL_API butter_ap(double rp, int ord, double *b, double *a)
Function calculates the transfer function coefficients of analog normalized lowpass Butterworth filt...
Definition: filter_ap.c:175
int DSPL_API ellip_ap(double rp, double rs, int ord, double *b, double *a)
Function calculates the transfer function coefficients of analog normalized lowpass elliptic filter ...
Definition: filter_ap.c:1538
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:497
#define ERROR_FILTER_ORD_BP
The order of the bandpass or notch filter is not set correctly. The order of the bandpass and notch f...
Definition: dspl.h:515