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 
100 #endif
101 #ifdef DOXYGEN_RUSSIAN
102 
174 #endif
175 int DSPL_API butter_ap(double rp, int ord, double* b, double* a)
176 {
177  int res;
178  complex_t *z = NULL;
179  complex_t *p = NULL;
180  int nz, np;
181 
182  if(rp < 0.0)
183  return ERROR_FILTER_RP;
184  if(ord < 1)
185  return ERROR_FILTER_ORD;
186  if(!a || !b)
187  return ERROR_PTR;
188 
189  z = (complex_t*) malloc(ord*sizeof(complex_t));
190  p = (complex_t*) malloc(ord*sizeof(complex_t));
191 
192 
193  res = butter_ap_zp(ord, rp, z, &nz, p, &np);
194  if(res != RES_OK)
195  goto exit_label;
196 
197  res = filter_zp2ab(z, nz, p, np, ord, b, a);
198  if(res != RES_OK)
199  goto exit_label;
200 
201  b[0] = a[0];
202 
203 
204 exit_label:
205  if(z)
206  free(z);
207  if(p)
208  free(p);
209  return res;
210 }
211 
212 
213 
214 
215 
216 
217 
218 
219 #ifdef DOXYGEN_ENGLISH
220 
305 #endif
306 #ifdef DOXYGEN_RUSSIAN
307 
401 #endif
402 int DSPL_API butter_ap_zp(int ord, double rp, complex_t* z, int* nz,
403  complex_t *p, int* np)
404 {
405  double alpha;
406  double theta;
407  double ep;
408  int r;
409  int L;
410  int ind = 0, k;
411 
412  if(rp < 0 || rp == 0)
413  return ERROR_FILTER_RP;
414  if(ord < 1)
415  return ERROR_FILTER_ORD;
416  if(!z || !p || !nz || !np)
417  return ERROR_PTR;
418 
419  ep = sqrt(pow(10.0, rp*0.1) - 1.0);
420  r = ord % 2;
421  L = (int)((ord-r)/2);
422 
423  alpha = pow(ep, -1.0/(double)ord);
424  if(r)
425  {
426  RE(p[ind]) = -alpha;
427  IM(p[ind]) = 0.0;
428  ind++;
429  }
430  for(k = 0; k < L; k++)
431  {
432  theta = M_PI*(double)(2*k + 1)/(double)(2*ord);
433  RE(p[ind]) = RE(p[ind+1]) = -alpha * sin(theta);
434  IM(p[ind]) = alpha * cos(theta);
435  IM(p[ind+1]) = -alpha * cos(theta);
436  ind+=2;
437  }
438  *np = ord;
439  *nz = 0;
440  return RES_OK;
441 }
442 
443 
444 
445 
446 
447 
448 #ifdef DOXYGEN_ENGLISH
449 
520 #endif
521 #ifdef DOXYGEN_RUSSIAN
522 
600 #endif
601 int DSPL_API cheby1_ap(double rp, int ord, double* b, double* a)
602 {
603  int res;
604  complex_t *z = NULL;
605  complex_t *p = NULL;
606  int nz, np, k;
607  complex_t h0 = {1.0, 0.0};
608  double tmp;
609 
610 
611  if(rp < 0.0)
612  return ERROR_FILTER_RP;
613  if(ord < 1)
614  return ERROR_FILTER_ORD;
615  if(!a || !b)
616  return ERROR_PTR;
617 
618  z = (complex_t*) malloc(ord*sizeof(complex_t));
619  p = (complex_t*) malloc(ord*sizeof(complex_t));
620 
621 
622  res = cheby1_ap_zp(ord, rp, z, &nz, p, &np);
623  if(res != RES_OK)
624  goto exit_label;
625 
626  res = filter_zp2ab(z, nz, p, np, ord, b, a);
627  if(res != RES_OK)
628  goto exit_label;
629 
630 
631  if(!(ord % 2))
632  RE(h0) = 1.0 / pow(10.0, rp*0.05);
633 
634  for(k = 0; k < np; k++)
635  {
636  tmp = CMRE(h0, p[k]);
637  IM(h0) = CMIM(h0, p[k]);
638  RE(h0) = tmp;
639  }
640 
641  b[0] = fabs(RE(h0));
642 
643 exit_label:
644  if(z)
645  free(z);
646  if(p)
647  free(p);
648  return res;
649 }
650 
651 
652 
653 
654 
655 #ifdef DOXYGEN_ENGLISH
656 
741 #endif
742 #ifdef DOXYGEN_RUSSIAN
743 
831 #endif
832 int DSPL_API cheby1_ap_zp(int ord, double rp, complex_t* z, int* nz,
833  complex_t* p, int* np)
834 {
835  double theta;
836  double ep;
837  double beta;
838  double shbeta;
839  double chbeta;
840  int r;
841  int L;
842  int ind = 0, k;
843 
844  if(rp < 0 || rp == 0)
845  return ERROR_FILTER_RP;
846  if(ord < 1)
847  return ERROR_FILTER_ORD;
848  if(!z || !p || !nz || !np)
849  return ERROR_PTR;
850 
851  ep = sqrt(pow(10.0, rp*0.1) - 1.0);
852  r = ord % 2;
853  L = (int)((ord-r)/2);
854 
855 
856  beta = asinh(1.0/ep)/(double)ord;
857  chbeta = cosh(beta);
858  shbeta = sinh(beta);
859 
860  if(r)
861  {
862  RE(p[ind]) = -shbeta;
863  IM(p[ind]) = 0.0;
864  ind++;
865  }
866  for(k = 0; k < L; k++)
867  {
868  theta = M_PI*(double)(2*k + 1)/(double)(2*ord);
869  RE(p[ind]) = RE(p[ind+1]) = -shbeta * sin(theta);
870  IM(p[ind]) = chbeta * cos(theta);
871  IM(p[ind+1]) = -IM(p[ind]);
872  ind+=2;
873  }
874  *np = ord;
875  *nz = 0;
876  return RES_OK;
877 }
878 
879 
880 
881 #ifdef DOXYGEN_ENGLISH
882 
955 #endif
956 #ifdef DOXYGEN_RUSSIAN
957 
1036 #endif
1037 int DSPL_API cheby2_ap(double rs, int ord, double* b, double* a)
1038 {
1039  int res;
1040  complex_t *z = NULL;
1041  complex_t *p = NULL;
1042  int nz, np;
1043  double norm;
1044 
1045 
1046  if(rs < 0.0)
1047  return ERROR_FILTER_RP;
1048  if(ord < 1)
1049  return ERROR_FILTER_ORD;
1050  if(!a || !b)
1051  return ERROR_PTR;
1052 
1053  z = (complex_t*) malloc(ord*sizeof(complex_t));
1054  p = (complex_t*) malloc(ord*sizeof(complex_t));
1055 
1056 
1057  res = cheby2_ap_zp(ord, rs, z, &nz, p, &np);
1058  if(res != RES_OK)
1059  goto exit_label;
1060 
1061  res = filter_zp2ab(z, nz, p, np, ord, b, a);
1062  if(res != RES_OK)
1063  goto exit_label;
1064 
1065  norm = a[0] / b[0];
1066 
1067  for(nz = 0; nz < ord+1; nz++)
1068  b[nz]*=norm;
1069 
1070 exit_label:
1071  if(z)
1072  free(z);
1073  if(p)
1074  free(p);
1075  return res;
1076 }
1077 
1078 
1079 
1080 #ifdef DOXYGEN_ENGLISH
1081 
1082 #endif
1083 #ifdef DOXYGEN_RUSSIAN
1084 
1085 #endif
1086 int DSPL_API cheby2_ap_wp1(double rp, double rs, int ord, double* b, double* a)
1087 {
1088  int err;
1089  double es, gp, alpha, beta, y, wp;
1090 
1091  if(rp <= 0)
1092  return ERROR_FILTER_RP;
1093 
1094  err = cheby2_ap(rs, ord, b, a);
1095  if(err!=RES_OK)
1096  goto exit_label;
1097 
1098  es = sqrt(pow(10.0, rs*0.1) - 1.0);
1099  gp = pow(10.0, -rp*0.05);
1100  alpha = gp * es / sqrt(1.0 - gp*gp);
1101  beta = alpha + sqrt(alpha * alpha - 1.0);
1102  y = log(beta)/ (double)ord;
1103  wp = 2.0 / (exp(y) + exp(-y));
1104 
1105  err = low2low(b, a, ord, wp, 1.0, b, a);
1106 
1107 exit_label:
1108  return err;
1109 }
1110 
1111 
1112 
1113 
1114 #ifdef DOXYGEN_ENGLISH
1115 
1206 #endif
1207 #ifdef DOXYGEN_RUSSIAN
1208 
1303 #endif
1304 int DSPL_API cheby2_ap_zp(int ord, double rs, complex_t* z, int* nz,
1305  complex_t *p, int* np)
1306 {
1307  double es;
1308  int L, r, k;
1309  double beta;
1310  int iz, ip;
1311 
1312  double alpha;
1313  double chb, shb, sa, ca;
1314  double ssh2, cch2;
1315 
1316  if(rs < 0 || rs == 0)
1317  return ERROR_FILTER_RS;
1318  if(ord < 1)
1319  return ERROR_FILTER_ORD;
1320  if(!z || !p || !nz || !np)
1321  return ERROR_PTR;
1322 
1323  es = sqrt(pow(10.0, rs*0.1) - 1.0);
1324  r = ord % 2;
1325  L = (int)((ord-r)/2);
1326 
1327  beta = asinh(es)/(double)ord;
1328 
1329  chb = cosh(beta);
1330  shb = sinh(beta);
1331 
1332  iz = ip = 0;
1333 
1334  if(r)
1335  {
1336  RE(p[0]) = -1.0 / sinh(beta);
1337  IM(p[0]) = 0.0;
1338  ip = 1;
1339  }
1340 
1341  for(k = 0; k < L; k++)
1342  {
1343  alpha = M_PI*(double)(2*k + 1)/(double)(2*ord);
1344  sa = sin(alpha);
1345  ca = cos(alpha);
1346  ssh2 = sa*shb;
1347  ssh2 *= ssh2;
1348 
1349  cch2 = ca*chb;
1350  cch2 *= cch2;
1351 
1352  RE(z[iz]) = RE(z[iz+1]) = 0.0;
1353  IM(z[iz]) = 1.0 / ca;
1354  IM(z[iz+1]) = -IM(z[iz]);
1355  iz+=2;
1356 
1357  RE(p[ip]) = RE(p[ip+1]) = -sa*shb / (ssh2 + cch2);
1358  IM(p[ip]) = ca*chb / (ssh2 + cch2);
1359  IM(p[ip+1]) = -IM(p[ip]);
1360  ip+=2;
1361  }
1362  *nz = iz;
1363  *np = ip;
1364 
1365  return RES_OK;
1366 }
1367 
1368 
1369 
1370 
1371 
1372 
1373 
1374 #ifdef DOXYGEN_ENGLISH
1375 
1451 #endif
1452 #ifdef DOXYGEN_RUSSIAN
1453 
1537 #endif
1538 int DSPL_API ellip_ap(double rp, double rs, int ord, double* b, double* a)
1539 {
1540  int res;
1541  complex_t *z = NULL;
1542  complex_t *p = NULL;
1543  int nz, np;
1544  double norm, g0;
1545 
1546 
1547  if(rp < 0.0)
1548  return ERROR_FILTER_RP;
1549  if(rs < 0.0)
1550  return ERROR_FILTER_RS;
1551  if(ord < 1)
1552  return ERROR_FILTER_ORD;
1553  if(!a || !b)
1554  return ERROR_PTR;
1555 
1556  z = (complex_t*) malloc(ord*sizeof(complex_t));
1557  p = (complex_t*) malloc(ord*sizeof(complex_t));
1558 
1559 
1560  res = ellip_ap_zp(ord, rp, rs, z, &nz, p, &np);
1561  if(res != RES_OK)
1562  goto exit_label;
1563 
1564  res = filter_zp2ab(z, nz, p, np, ord, b, a);
1565  if(res != RES_OK)
1566  goto exit_label;
1567 
1568 
1569  g0 = 1.0;
1570  if(!(ord % 2))
1571  {
1572  g0 = 1.0 / pow(10.0, rp*0.05);
1573  }
1574 
1575 
1576  norm = g0 * a[0] / b[0];
1577 
1578  for(nz = 0; nz < ord+1; nz++)
1579  b[nz]*=norm;
1580 
1581  exit_label:
1582  if(z)
1583  free(z);
1584  if(p)
1585  free(p);
1586  return res;
1587 }
1588 
1589 
1590 
1591 
1592 
1593 
1594 #ifdef DOXYGEN_ENGLISH
1595 
1686 #endif
1687 #ifdef DOXYGEN_RUSSIAN
1688 
1785 #endif
1786 int DSPL_API ellip_ap_zp(int ord, double rp, double rs,
1787  complex_t* z, int* nz, complex_t* p, int* np)
1788 {
1789  double es, ep;
1790  int L, r, n, res;
1791  int iz, ip;
1792  double ke, k, u, t;
1793  complex_t tc, v0, jv0;
1794 
1795 
1796  if(rp < 0 || rp == 0)
1797  return ERROR_FILTER_RP;
1798  if(rs < 0 || rs == 0)
1799  return ERROR_FILTER_RS;
1800  if(ord < 1)
1801  return ERROR_FILTER_ORD;
1802  if(!z || !p || !nz || !np)
1803  return ERROR_PTR;
1804 
1805  es = sqrt(pow(10.0, rs*0.1) - 1.0);
1806  ep = sqrt(pow(10.0, rp*0.1) - 1.0);
1807  ke = ep / es;
1808 
1809  r = ord % 2;
1810  L = (int)((ord-r)/2);
1811 
1812  res = ellip_modulareq(rp, rs, ord, &k);
1813  if(res != RES_OK)
1814  return res;
1815  // v0
1816  RE(tc) = 0.0;
1817  IM(tc) = 1.0 / ep;
1818 
1819  ellip_asn_cmplx(&tc, 1, ke, &v0);
1820 
1821  t = RE(v0);
1822  RE(v0) = IM(v0) / (double)ord;
1823  IM(v0) = -t / (double)ord;
1824 
1825  RE(jv0) = -IM(v0);
1826  IM(jv0) = RE(v0);
1827 
1828  iz = ip = 0;
1829 
1830  if(r)
1831  {
1832  res = ellip_sn_cmplx(&jv0, 1, k, &tc);
1833  if(res != RES_OK)
1834  return res;
1835  RE(p[0]) = -IM(tc);
1836  IM(p[0]) = RE(tc);
1837  ip = 1;
1838  }
1839 
1840  for(n = 0; n < L; n++)
1841  {
1842  u = (double)(2 * n + 1)/(double)ord;
1843 
1844  res = ellip_cd(& u, 1, k, &t);
1845  if(res != RES_OK)
1846  return res;
1847 
1848  RE(z[iz]) = RE(z[iz+1]) = 0.0;
1849  IM(z[iz]) = 1.0/(k*t);
1850  IM(z[iz+1]) = -1.0/(k*t);
1851  iz+=2;
1852 
1853  RE(tc) = u - RE(jv0);
1854  IM(tc) = - IM(jv0);
1855 
1856  res = ellip_cd_cmplx(&tc, 1, k, p+ip+1);
1857  if(res != RES_OK)
1858  return res;
1859 
1860  RE(p[ip]) = -IM(p[ip+1]);
1861  IM(p[ip]) = RE(p[ip+1]);
1862 
1863  RE(p[ip+1]) = RE(p[ip]);
1864  IM(p[ip+1]) = -IM(p[ip]);
1865 
1866  ip+=2;
1867  }
1868  *nz = iz;
1869  *np = ip;
1870 
1871  return RES_OK;
1872 }
1873 
1874 
1875 
1876 
1877 #ifdef DOXYGEN_ENGLISH
1878 
1946 #endif
1947 #ifdef DOXYGEN_RUSSIAN
1948 
2015 #endif
2016 int DSPL_API filter_zp2ab(complex_t* z, int nz, complex_t* p, int np,
2017  int ord, double* b, double* a)
2018 {
2019  complex_t *acc = NULL;
2020  int res;
2021 
2022  if(!z || !p || !b || !a)
2023  return ERROR_PTR;
2024  if(nz < 0 || np < 0)
2025  return ERROR_SIZE;
2026  if(nz > ord || np > ord)
2027  return ERROR_POLY_ORD;
2028 
2029  acc = (complex_t*) malloc((ord+1) * sizeof(complex_t));
2030  res = poly_z2a_cmplx(z, nz, ord, acc);
2031  if(res != RES_OK)
2032  goto exit_label;
2033 
2034  res = cmplx2re(acc, ord+1, b, NULL);
2035  if(res != RES_OK)
2036  goto exit_label;
2037 
2038  res = poly_z2a_cmplx(p, np, ord, acc);
2039  if(res != RES_OK)
2040  goto exit_label;
2041 
2042  res = cmplx2re(acc, ord+1, a, NULL);
2043  if(res != RES_OK)
2044  goto exit_label;
2045 
2046 exit_label:
2047  if(acc)
2048  free(acc);
2049  return res;
2050 
2051 }
2052 
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:1786
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 ellip_cd(double *u, int n, double k, double *y)
Jacobi elliptic function of real vector argument.
Definition: ellipj.c:577
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:402
#define ERROR_POLY_ORD
The polynomial order is incorrect. The order of the polynomial must be a positive integer.
Definition: dspl.h:544
int DSPL_API low2low(double *b, double *a, int ord, double w0, double w1, double *beta, double *alpha)
Definition: filter_ft.c:426
#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:514
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:359
#define ERROR_PTR
Pointer error. This error means that one of the required pointers (memory to be allocated for) is tra...
Definition: dspl.h:545
#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:553
#define ERROR_FILTER_RP
The filter non-uniformity parameter in the passband is set incorrectly. This parameter is specified i...
Definition: dspl.h:516
#define ERROR_FILTER_RS
The filter suppression parameter in the boom bar is not set correctly. This parameter is specified in...
Definition: dspl.h:517
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:1179
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
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 cmplx2re(complex_t *x, int n, double *re, double *im)
Separate complex vector to the real and image vectors.
Definition: complex.c:130
double complex_t[2]
Complex data type.
Definition: dspl.h:86
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:832
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
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:684
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:417
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:1304
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:1037