libdspl-2.0
Digital Signal Processing Algorithm Library
filter_ap.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 <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include "dspl.h"
25 
26 
27 
28 #ifdef DOXYGEN_ENGLISH
29 
99 #endif
100 #ifdef DOXYGEN_RUSSIAN
101 
173 #endif
174 int DSPL_API butter_ap(double rp, int ord, double* b, double* a)
175 {
176  int res;
177  complex_t *z = NULL;
178  complex_t *p = NULL;
179  int nz, np;
180 
181  if(rp < 0.0)
182  return ERROR_FILTER_RP;
183  if(ord < 1)
184  return ERROR_FILTER_ORD;
185  if(!a || !b)
186  return ERROR_PTR;
187 
188  z = (complex_t*) malloc(ord*sizeof(complex_t));
189  p = (complex_t*) malloc(ord*sizeof(complex_t));
190 
191 
192  res = butter_ap_zp(ord, rp, z, &nz, p, &np);
193  if(res != RES_OK)
194  goto exit_label;
195 
196  res = filter_zp2ab(z, nz, p, np, ord, b, a);
197  if(res != RES_OK)
198  goto exit_label;
199 
200  b[0] = a[0];
201 
202 
203 exit_label:
204  if(z)
205  free(z);
206  if(p)
207  free(p);
208  return res;
209 }
210 
211 
212 
213 #ifdef DOXYGEN_ENGLISH
214 
299 #endif
300 #ifdef DOXYGEN_RUSSIAN
301 
395 #endif
396 int DSPL_API butter_ap_zp(int ord, double rp, complex_t* z, int* nz,
397  complex_t *p, int* np)
398 {
399  double alpha;
400  double theta;
401  double ep;
402  int r;
403  int L;
404  int ind = 0, k;
405 
406  if(rp < 0 || rp == 0)
407  return ERROR_FILTER_RP;
408  if(ord < 1)
409  return ERROR_FILTER_ORD;
410  if(!z || !p || !nz || !np)
411  return ERROR_PTR;
412 
413  ep = sqrt(pow(10.0, rp*0.1) - 1.0);
414  r = ord % 2;
415  L = (int)((ord-r)/2);
416 
417  alpha = pow(ep, -1.0/(double)ord);
418  if(r)
419  {
420  RE(p[ind]) = -alpha;
421  IM(p[ind]) = 0.0;
422  ind++;
423  }
424  for(k = 0; k < L; k++)
425  {
426  theta = M_PI*(double)(2*k + 1)/(double)(2*ord);
427  RE(p[ind]) = RE(p[ind+1]) = -alpha * sin(theta);
428  IM(p[ind]) = alpha * cos(theta);
429  IM(p[ind+1]) = -alpha * cos(theta);
430  ind+=2;
431  }
432  *np = ord;
433  *nz = 0;
434  return RES_OK;
435 }
436 
437 
438 
439 
440 
441 
442 #ifdef DOXYGEN_ENGLISH
443 
514 #endif
515 #ifdef DOXYGEN_RUSSIAN
516 
594 #endif
595 int DSPL_API cheby1_ap(double rp, int ord, double* b, double* a)
596 {
597  int res;
598  complex_t *z = NULL;
599  complex_t *p = NULL;
600  int nz, np, k;
601  complex_t h0 = {1.0, 0.0};
602  double tmp;
603 
604 
605  if(rp < 0.0)
606  return ERROR_FILTER_RP;
607  if(ord < 1)
608  return ERROR_FILTER_ORD;
609  if(!a || !b)
610  return ERROR_PTR;
611 
612  z = (complex_t*) malloc(ord*sizeof(complex_t));
613  p = (complex_t*) malloc(ord*sizeof(complex_t));
614 
615 
616  res = cheby1_ap_zp(ord, rp, z, &nz, p, &np);
617  if(res != RES_OK)
618  goto exit_label;
619 
620  res = filter_zp2ab(z, nz, p, np, ord, b, a);
621  if(res != RES_OK)
622  goto exit_label;
623 
624 
625  if(!(ord % 2))
626  RE(h0) = 1.0 / pow(10.0, rp*0.05);
627 
628  for(k = 0; k < np; k++)
629  {
630  tmp = CMRE(h0, p[k]);
631  IM(h0) = CMIM(h0, p[k]);
632  RE(h0) = tmp;
633  }
634 
635  b[0] = fabs(RE(h0));
636 
637 exit_label:
638  if(z)
639  free(z);
640  if(p)
641  free(p);
642  return res;
643 }
644 
645 
646 
647 
648 
649 #ifdef DOXYGEN_ENGLISH
650 
735 #endif
736 #ifdef DOXYGEN_RUSSIAN
737 
825 #endif
826 int DSPL_API cheby1_ap_zp(int ord, double rp, complex_t* z, int* nz,
827  complex_t* p, int* np)
828 {
829  double theta;
830  double ep;
831  double beta;
832  double shbeta;
833  double chbeta;
834  int r;
835  int L;
836  int ind = 0, k;
837 
838  if(rp < 0 || rp == 0)
839  return ERROR_FILTER_RP;
840  if(ord < 1)
841  return ERROR_FILTER_ORD;
842  if(!z || !p || !nz || !np)
843  return ERROR_PTR;
844 
845  ep = sqrt(pow(10.0, rp*0.1) - 1.0);
846  r = ord % 2;
847  L = (int)((ord-r)/2);
848 
849 
850  beta = asinh(1.0/ep)/(double)ord;
851  chbeta = cosh(beta);
852  shbeta = sinh(beta);
853 
854  if(r)
855  {
856  RE(p[ind]) = -shbeta;
857  IM(p[ind]) = 0.0;
858  ind++;
859  }
860  for(k = 0; k < L; k++)
861  {
862  theta = M_PI*(double)(2*k + 1)/(double)(2*ord);
863  RE(p[ind]) = RE(p[ind+1]) = -shbeta * sin(theta);
864  IM(p[ind]) = chbeta * cos(theta);
865  IM(p[ind+1]) = -IM(p[ind]);
866  ind+=2;
867  }
868  *np = ord;
869  *nz = 0;
870  return RES_OK;
871 }
872 
873 
874 
875 #ifdef DOXYGEN_ENGLISH
876 
949 #endif
950 #ifdef DOXYGEN_RUSSIAN
951 
1030 #endif
1031 int DSPL_API cheby2_ap(double rs, int ord, double* b, double* a)
1032 {
1033  int res;
1034  complex_t *z = NULL;
1035  complex_t *p = NULL;
1036  int nz, np;
1037  double norm;
1038 
1039 
1040  if(rs < 0.0)
1041  return ERROR_FILTER_RP;
1042  if(ord < 1)
1043  return ERROR_FILTER_ORD;
1044  if(!a || !b)
1045  return ERROR_PTR;
1046 
1047  z = (complex_t*) malloc(ord*sizeof(complex_t));
1048  p = (complex_t*) malloc(ord*sizeof(complex_t));
1049 
1050 
1051  res = cheby2_ap_zp(ord, rs, z, &nz, p, &np);
1052  if(res != RES_OK)
1053  goto exit_label;
1054 
1055  res = filter_zp2ab(z, nz, p, np, ord, b, a);
1056  if(res != RES_OK)
1057  goto exit_label;
1058 
1059  norm = a[0] / b[0];
1060 
1061  for(nz = 0; nz < ord+1; nz++)
1062  b[nz]*=norm;
1063 
1064 exit_label:
1065  if(z)
1066  free(z);
1067  if(p)
1068  free(p);
1069  return res;
1070 }
1071 
1072 
1073 
1074 #ifdef DOXYGEN_ENGLISH
1075 
1076 #endif
1077 #ifdef DOXYGEN_RUSSIAN
1078 
1079 #endif
1080 int DSPL_API cheby2_ap_wp1(double rp, double rs, int ord, double* b, double* a)
1081 {
1082  int err;
1083  double es, gp, alpha, beta, y, wp;
1084 
1085  if(rp <= 0)
1086  return ERROR_FILTER_RP;
1087 
1088  err = cheby2_ap(rs, ord, b, a);
1089  if(err!=RES_OK)
1090  goto exit_label;
1091 
1092  es = sqrt(pow(10.0, rs*0.1) - 1.0);
1093  gp = pow(10.0, -rp*0.05);
1094  alpha = gp * es / sqrt(1.0 - gp*gp);
1095  beta = alpha + sqrt(alpha * alpha - 1.0);
1096  y = log(beta)/ (double)ord;
1097  wp = 2.0 / (exp(y) + exp(-y));
1098 
1099  err = low2low(b, a, ord, wp, 1.0, b, a);
1100 
1101 exit_label:
1102  return err;
1103 }
1104 
1105 
1106 
1107 
1108 #ifdef DOXYGEN_ENGLISH
1109 
1200 #endif
1201 #ifdef DOXYGEN_RUSSIAN
1202 
1297 #endif
1298 int DSPL_API cheby2_ap_zp(int ord, double rs, complex_t* z, int* nz,
1299  complex_t *p, int* np)
1300 {
1301  double es;
1302  int L, r, k;
1303  double beta;
1304  int iz, ip;
1305 
1306  double alpha;
1307  double chb, shb, sa, ca;
1308  double ssh2, cch2;
1309 
1310  if(rs < 0 || rs == 0)
1311  return ERROR_FILTER_RS;
1312  if(ord < 1)
1313  return ERROR_FILTER_ORD;
1314  if(!z || !p || !nz || !np)
1315  return ERROR_PTR;
1316 
1317  es = sqrt(pow(10.0, rs*0.1) - 1.0);
1318  r = ord % 2;
1319  L = (int)((ord-r)/2);
1320 
1321  beta = asinh(es)/(double)ord;
1322 
1323  chb = cosh(beta);
1324  shb = sinh(beta);
1325 
1326  iz = ip = 0;
1327 
1328  if(r)
1329  {
1330  RE(p[0]) = -1.0 / sinh(beta);
1331  IM(p[0]) = 0.0;
1332  ip = 1;
1333  }
1334 
1335  for(k = 0; k < L; k++)
1336  {
1337  alpha = M_PI*(double)(2*k + 1)/(double)(2*ord);
1338  sa = sin(alpha);
1339  ca = cos(alpha);
1340  ssh2 = sa*shb;
1341  ssh2 *= ssh2;
1342 
1343  cch2 = ca*chb;
1344  cch2 *= cch2;
1345 
1346  RE(z[iz]) = RE(z[iz+1]) = 0.0;
1347  IM(z[iz]) = 1.0 / ca;
1348  IM(z[iz+1]) = -IM(z[iz]);
1349  iz+=2;
1350 
1351  RE(p[ip]) = RE(p[ip+1]) = -sa*shb / (ssh2 + cch2);
1352  IM(p[ip]) = ca*chb / (ssh2 + cch2);
1353  IM(p[ip+1]) = -IM(p[ip]);
1354  ip+=2;
1355  }
1356  *nz = iz;
1357  *np = ip;
1358 
1359  return RES_OK;
1360 }
1361 
1362 
1363 
1364 
1365 
1366 
1367 
1368 #ifdef DOXYGEN_ENGLISH
1369 
1445 #endif
1446 #ifdef DOXYGEN_RUSSIAN
1447 
1531 #endif
1532 int DSPL_API ellip_ap(double rp, double rs, int ord, double* b, double* a)
1533 {
1534  int res;
1535  complex_t *z = NULL;
1536  complex_t *p = NULL;
1537  int nz, np;
1538  double norm, g0;
1539 
1540 
1541  if(rp < 0.0)
1542  return ERROR_FILTER_RP;
1543  if(rs < 0.0)
1544  return ERROR_FILTER_RS;
1545  if(ord < 1)
1546  return ERROR_FILTER_ORD;
1547  if(!a || !b)
1548  return ERROR_PTR;
1549 
1550  z = (complex_t*) malloc(ord*sizeof(complex_t));
1551  p = (complex_t*) malloc(ord*sizeof(complex_t));
1552 
1553 
1554  res = ellip_ap_zp(ord, rp, rs, z, &nz, p, &np);
1555  if(res != RES_OK)
1556  goto exit_label;
1557 
1558  res = filter_zp2ab(z, nz, p, np, ord, b, a);
1559  if(res != RES_OK)
1560  goto exit_label;
1561 
1562 
1563  g0 = 1.0;
1564  if(!(ord % 2))
1565  {
1566  g0 = 1.0 / pow(10.0, rp*0.05);
1567  }
1568 
1569 
1570  norm = g0 * a[0] / b[0];
1571 
1572  for(nz = 0; nz < ord+1; nz++)
1573  b[nz]*=norm;
1574 
1575  exit_label:
1576  if(z)
1577  free(z);
1578  if(p)
1579  free(p);
1580  return res;
1581 }
1582 
1583 
1584 
1585 
1586 
1587 
1588 #ifdef DOXYGEN_ENGLISH
1589 
1680 #endif
1681 #ifdef DOXYGEN_RUSSIAN
1682 
1779 #endif
1780 int DSPL_API ellip_ap_zp(int ord, double rp, double rs,
1781  complex_t* z, int* nz, complex_t* p, int* np)
1782 {
1783  double es, ep;
1784  int L, r, n, res;
1785  int iz, ip;
1786  double ke, k, u, t;
1787  complex_t tc, v0, jv0;
1788 
1789 
1790  if(rp < 0 || rp == 0)
1791  return ERROR_FILTER_RP;
1792  if(rs < 0 || rs == 0)
1793  return ERROR_FILTER_RS;
1794  if(ord < 1)
1795  return ERROR_FILTER_ORD;
1796  if(!z || !p || !nz || !np)
1797  return ERROR_PTR;
1798 
1799  es = sqrt(pow(10.0, rs*0.1) - 1.0);
1800  ep = sqrt(pow(10.0, rp*0.1) - 1.0);
1801  ke = ep / es;
1802 
1803  r = ord % 2;
1804  L = (int)((ord-r)/2);
1805 
1806  res = ellip_modulareq(rp, rs, ord, &k);
1807  if(res != RES_OK)
1808  return res;
1809  // v0
1810  RE(tc) = 0.0;
1811  IM(tc) = 1.0 / ep;
1812 
1813  ellip_asn_cmplx(&tc, 1, ke, &v0);
1814 
1815  t = RE(v0);
1816  RE(v0) = IM(v0) / (double)ord;
1817  IM(v0) = -t / (double)ord;
1818 
1819  RE(jv0) = -IM(v0);
1820  IM(jv0) = RE(v0);
1821 
1822  iz = ip = 0;
1823 
1824  if(r)
1825  {
1826  res = ellip_sn_cmplx(&jv0, 1, k, &tc);
1827  if(res != RES_OK)
1828  return res;
1829  RE(p[0]) = -IM(tc);
1830  IM(p[0]) = RE(tc);
1831  ip = 1;
1832  }
1833 
1834  for(n = 0; n < L; n++)
1835  {
1836  u = (double)(2 * n + 1)/(double)ord;
1837 
1838  res = ellip_cd(& u, 1, k, &t);
1839  if(res != RES_OK)
1840  return res;
1841 
1842  RE(z[iz]) = RE(z[iz+1]) = 0.0;
1843  IM(z[iz]) = 1.0/(k*t);
1844  IM(z[iz+1]) = -1.0/(k*t);
1845  iz+=2;
1846 
1847  RE(tc) = u - RE(jv0);
1848  IM(tc) = - IM(jv0);
1849 
1850  res = ellip_cd_cmplx(&tc, 1, k, p+ip+1);
1851  if(res != RES_OK)
1852  return res;
1853 
1854  RE(p[ip]) = -IM(p[ip+1]);
1855  IM(p[ip]) = RE(p[ip+1]);
1856 
1857  RE(p[ip+1]) = RE(p[ip]);
1858  IM(p[ip+1]) = -IM(p[ip]);
1859 
1860  ip+=2;
1861  }
1862  *nz = iz;
1863  *np = ip;
1864 
1865  return RES_OK;
1866 }
1867 
1868 
1869 
1870 
1871 #ifdef DOXYGEN_ENGLISH
1872 
1940 #endif
1941 #ifdef DOXYGEN_RUSSIAN
1942 
2009 #endif
2010 int DSPL_API filter_zp2ab(complex_t* z, int nz, complex_t* p, int np,
2011  int ord, double* b, double* a)
2012 {
2013  complex_t *acc = NULL;
2014  int res;
2015 
2016  if(!z || !p || !b || !a)
2017  return ERROR_PTR;
2018  if(nz < 0 || np < 0)
2019  return ERROR_SIZE;
2020  if(nz > ord || np > ord)
2021  return ERROR_POLY_ORD;
2022 
2023  acc = (complex_t*) malloc((ord+1) * sizeof(complex_t));
2024  res = poly_z2a_cmplx(z, nz, ord, acc);
2025  if(res != RES_OK)
2026  goto exit_label;
2027 
2028  res = cmplx2re(acc, ord+1, b, NULL);
2029  if(res != RES_OK)
2030  goto exit_label;
2031 
2032  res = poly_z2a_cmplx(p, np, ord, acc);
2033  if(res != RES_OK)
2034  goto exit_label;
2035 
2036  res = cmplx2re(acc, ord+1, a, NULL);
2037  if(res != RES_OK)
2038  goto exit_label;
2039 
2040 exit_label:
2041  if(acc)
2042  free(acc);
2043  return res;
2044 
2045 }
2046 
ellip_ap_zp
int DSPL_API ellip_ap_zp(int ord, double rp, double rs, complex_t *z, int *nz, complex_t *p, int *np)
Function calculates arrays of zeros and poles for analog normlized lowpass elliptic filter transfer f...
Definition: filter_ap.c:1780
cheby1_ap
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:595
ellip_cd
int DSPL_API ellip_cd(double *u, int n, double k, double *y)
Jacobi elliptic function of real vector argument.
Definition: ellipj.c:594
butter_ap_zp
int DSPL_API butter_ap_zp(int ord, double rp, complex_t *z, int *nz, complex_t *p, int *np)
Function calculates arrays of zeros and poles for analog normlized lowpass Batterworth filter transfe...
Definition: filter_ap.c:396
ERROR_POLY_ORD
#define ERROR_POLY_ORD
The polynomial order is incorrect. The order of the polynomial must be a positive integer.
Definition: dspl.h:609
low2low
int DSPL_API low2low(double *b, double *a, int ord, double w0, double w1, double *beta, double *alpha)
Definition: filter_ft.c:426
ERROR_FILTER_ORD
#define ERROR_FILTER_ORD
The filter order is incorrect. The order of the filter must be specified by a positive integer value.
Definition: dspl.h:575
RE
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:420
ERROR_PTR
#define ERROR_PTR
Pointer error. This error means that one of the required pointers (memory to be allocated for) is tra...
Definition: dspl.h:610
ERROR_SIZE
#define ERROR_SIZE
Error array size. This error occurs when in addition to the pointer the wrong input is passed to the ...
Definition: dspl.h:618
ERROR_FILTER_RP
#define ERROR_FILTER_RP
The filter non-uniformity parameter in the passband is set incorrectly. This parameter is specified i...
Definition: dspl.h:577
ERROR_FILTER_RS
#define ERROR_FILTER_RS
The filter suppression parameter in the boom bar is not set correctly. This parameter is specified in...
Definition: dspl.h:578
ellip_sn_cmplx
int DSPL_API ellip_sn_cmplx(complex_t *u, int n, double k, complex_t *y)
Jacobi elliptic function of complex vector argument.
Definition: ellipj.c:1216
ellip_asn_cmplx
int DSPL_API ellip_asn_cmplx(complex_t *w, int n, double k, complex_t *u)
Inverse Jacobi elliptic function of complex vector argument.
Definition: ellipj.c:447
butter_ap
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:174
cmplx2re
int DSPL_API cmplx2re(complex_t *x, int n, double *re, double *im)
Separate complex vector to the real and image vectors.
Definition: complex.c:130
complex_t
double complex_t[2]
Complex data type.
Definition: dspl.h:86
cheby1_ap_zp
int DSPL_API cheby1_ap_zp(int ord, double rp, complex_t *z, int *nz, complex_t *p, int *np)
Function calculates arrays of zeros and poles for analog normlized lowpass Chebyshev type 1 filter tr...
Definition: filter_ap.c:826
ellip_ap
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:1532
RES_OK
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:558
ellip_cd_cmplx
int DSPL_API ellip_cd_cmplx(complex_t *u, int n, double k, complex_t *y)
Jacobi elliptic function of complex vector argument.
Definition: ellipj.c:701
IM
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:478
cheby2_ap_zp
int DSPL_API cheby2_ap_zp(int ord, double rs, complex_t *z, int *nz, complex_t *p, int *np)
Function calculates arrays of zeros and poles for analog normlized lowpass Chebyshev type 2 filter tr...
Definition: filter_ap.c:1298
cheby2_ap
int DSPL_API cheby2_ap(double rs, int ord, double *b, double *a)
Function calculates the transfer function coefficients of analog normalized lowpass Chebyshev type 2...
Definition: filter_ap.c:1031
Generated on Mon Dec 14 2020 14:02:35 for libdspl-2.0 by   doxygen 1.8.20
libdspl-2.0 documentation is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.