libdspl-2.0
Digital Signal Processing Algorithm Library
filter_ft.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 
30 #endif
31 #ifdef DOXYGEN_RUSSIAN
32 
33 #endif
34 double DSPL_API filter_ws1(int ord, double rp, double rs, int type)
35 {
36  double es2, ep2, gs2, x, ws;
37 
38  if(ord<1 || rp < 0.0 || rs < 0.0)
39  return -1.0;
40 
41  es2 = pow(10.0, rs*0.1) - 1.0;
42  ep2 = pow(10.0, rp*0.1) - 1.0;
43  gs2 = 1.0 / (1.0 + es2);
44 
45  x = (1.0 - gs2) / (gs2 * ep2);
46 
47  switch( type & DSPL_FILTER_APPROX_MASK)
48  {
49  case DSPL_FILTER_BUTTER:
50  ws = pow(x, 0.5 / (double)ord);
51  break;
52  case DSPL_FILTER_CHEBY1:
53  case DSPL_FILTER_CHEBY2:
54  x = sqrt(x) + sqrt(x - 1.0);
55  x = log(x) / (double)ord;
56  ws = 0.5 * (exp(-x) + exp(x));
57  break;
58  case DSPL_FILTER_ELLIP:
59  {
60  double k, k1;
61  complex_t y, z;
62  int res;
63  k = sqrt(ep2 / es2);
64  res = ellip_modulareq(rp, rs, ord, &k1);
65  if(res != RES_OK)
66  {
67  ws = -1.0;
68  break;
69  }
70  RE(z) = sqrt(x);
71  IM(z) = 0.0;
72 
73  res = ellip_acd_cmplx(&z, 1, k, &y);
74  if(res != RES_OK)
75  {
76  ws = -1.0;
77  break;
78  }
79  RE(y) /= (double)ord;
80  IM(y) /= (double)ord;
81  res = ellip_cd_cmplx(&y, 1, k1, &z);
82  if(res != RES_OK)
83  {
84  ws = -1.0;
85  break;
86  }
87  ws = RE(z);
88  break;
89  }
90  default:
91  ws = -1.0;
92  break;
93  }
94  return ws;
95 }
96 
97 
98 #ifdef DOXYGEN_ENGLISH
99 
100 #endif
101 #ifdef DOXYGEN_RUSSIAN
102 
103 #endif
104 int DSPL_API low2bp(double* b, double* a, int ord,
105  double w0, double wpl, double wph,
106  double* beta, double* alpha)
107 {
108 
109  double num[3] = {0.0, 0.0, 1.0};
110  double den[3] = {0.0, 0.0, 0.0};
111 
112  if(!b || !a || !beta || !alpha)
113  return ERROR_PTR;
114  if(ord < 1)
115  return ERROR_FILTER_ORD;
116  if(w0 <= 0.0 || wpl <= 0.0 || wph <= 0.0 || wph <= wpl)
117  return ERROR_FILTER_FT;
118 
119  num[0] = (wph * wpl) / (w0 * w0);
120  den[1] = (wph - wpl) / w0;
121 
122  return ratcompos(b, a, ord, num, den, 2, beta, alpha);
123 }
124 
125 
126 
127 
128 
129 
130 #ifdef DOXYGEN_ENGLISH
131 
132 #endif
133 #ifdef DOXYGEN_RUSSIAN
134 
135 #endif
136 int DSPL_API low2bs(double* b, double* a, int ord,
137  double w0, double wsl, double wsh,
138  double* beta, double* alpha)
139 {
140 
141  double den[3] = {0.0, 0.0, 1.0};
142  double num[3] = {0.0, 0.0, 0.0};
143 
144  if(!b || !a || !beta || !alpha)
145  return ERROR_PTR;
146  if(ord < 1)
147  return ERROR_FILTER_ORD;
148  if(w0 <= 0.0 || wsl <= 0.0 || wsh <= 0.0 || wsh <= wsl)
149  return ERROR_FILTER_FT;
150 
151  den[0] = (wsh * wsl) / (w0 * w0);
152  num[1] = (wsh - wsl) / w0;
153 
154  return ratcompos(b, a, ord, num, den, 2, beta, alpha);
155 }
156 
157 
158 
159 
160 
161 #ifdef DOXYGEN_ENGLISH
162 
218 #endif
219 #ifdef DOXYGEN_RUSSIAN
220 
279 #endif
280 int DSPL_API low2high(double* b, double* a, int ord, double w0, double w1,
281  double* beta, double* alpha)
282 {
283 
284  double num[2] = {0.0, 0.0};
285  double den[2] = {0.0, 1.0};
286 
287  if(!b || !a || !beta || !alpha)
288  return ERROR_PTR;
289  if(ord < 1)
290  return ERROR_FILTER_ORD;
291  if(w0 <= 0.0 || w1 <= 0.0)
292  return ERROR_FILTER_FT;
293 
294  num[0] = w1 / w0;
295 
296  return ratcompos(b, a, ord, num, den, 1, beta, alpha);
297 }
298 
299 
300 
301 
302 
303 
304 
305 
306 #ifdef DOXYGEN_ENGLISH
307 
365 #endif
366 #ifdef DOXYGEN_RUSSIAN
367 
425 #endif
426 int DSPL_API low2low(double* b, double* a, int ord, double w0, double w1,
427  double* beta, double* alpha)
428 {
429 
430  double num[2] = {0.0, 1.0};
431  double den[2] = {0.0, 0.0};
432 
433  if(!b || !a || !beta || !alpha)
434  return ERROR_PTR;
435  if(ord < 1)
436  return ERROR_FILTER_ORD;
437  if(w0 <= 0.0 || w1 <= 0.0)
438  return ERROR_FILTER_FT;
439 
440  den[0] = w1 / w0;
441 
442  return ratcompos(b, a, ord, num, den, 1, beta, alpha);
443 }
444 
445 
446 
447 
448 #ifdef DOXYGEN_ENGLISH
449 
522 #endif
523 #ifdef DOXYGEN_RUSSIAN
524 
601 #endif
602 int DSPL_API ratcompos(double* b, double* a, int n,
603  double* c, double* d, int p,
604  double* beta, double* alpha)
605 {
606 
607  int k2, i, k, pn, pd, ln, ld, k2s, nk2s;
608  double *num = NULL, *den = NULL, *ndn = NULL, *ndd = NULL;
609  int res;
610 
611  if (!a || !b || !c || !d || !beta || !alpha)
612  {
613  res = ERROR_PTR;
614  goto exit_label;
615  }
616  if(n < 1 || p < 1)
617  {
618  res = ERROR_SIZE;
619  goto exit_label;
620  }
621 
622  k2 = (n*p)+1;
623  k2s = k2*sizeof(double); /* alpha and beta size */
624  nk2s = (n+1)*k2*sizeof(double); /* num, den, ndn and ndd size */
625 
626  num = (double*)malloc(nk2s);
627  den = (double*)malloc(nk2s);
628  ndn = (double*)malloc(nk2s);
629  ndd = (double*)malloc(nk2s);
630 
631  memset(num, 0, nk2s);
632  memset(den, 0, nk2s);
633  memset(ndn, 0, nk2s);
634  memset(ndd, 0, nk2s);
635 
636 
637  num[0] = den[0] = 1.0;
638  pn = 0;
639  ln = 1;
640  for(i = 1; i < n+1; i++)
641  {
642  res = conv(num+pn, ln, c, p+1, num+pn+k2);
643  if(res!=RES_OK)
644  goto exit_label;
645  res = conv(den+pn, ln, d, p+1, den+pn+k2);
646  if(res!=RES_OK)
647  goto exit_label;
648  pn += k2;
649  ln += p;
650  }
651 
652  pn = 0;
653  pd = n*k2;
654  ln = 1;
655  ld = k2;
656 
657  for (i = 0; i < n+1; i++)
658  {
659  res = conv(num + pn, ln, den + pd, ld, ndn + i*k2);
660  if(res!=RES_OK)
661  goto exit_label;
662  ln += p;
663  ld -= p;
664  pn += k2;
665  pd -= k2;
666  }
667 
668  for (i = 0; i < n+1; i++)
669  {
670  for (k = 0; k < k2; k++)
671  {
672  ndd[i*k2 + k] = ndn[i*k2 + k] * a[i];
673  ndn[i*k2 + k] *= b[i];
674  }
675  }
676 
677 
678  memset(alpha, 0, k2s);
679  memset(beta, 0, k2s);
680 
681  for (k = 0; k < k2; k++)
682  {
683  for (i = 0; i < n+1; i++)
684  {
685  beta[k] += ndn[i*k2 + k];
686  alpha[k] += ndd[i*k2 + k];
687  }
688  }
689 
690  res = RES_OK;
691 
692 exit_label:
693  if(num)
694  free(num);
695  if(den)
696  free(den);
697  if(ndn)
698  free(ndn);
699  if(ndd)
700  free(ndd);
701 
702  return res;
703 }
704 
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 conv(double *a, int na, double *b, int nb, double *c)
Real vectors linear convolution.
Definition: conv.c:157
#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
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_FT
The conversion frequencies of the low-pass filter and low-pass filter are incorrect....
Definition: dspl.h:513
double complex_t[2]
Complex data type.
Definition: dspl.h:86
#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
int DSPL_API ellip_acd_cmplx(complex_t *w, int n, double k, complex_t *u)
Inverse Jacobi elliptic function of complex vector argument.
Definition: ellipj.c:211
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:417