libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
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  * Recalculate ws frequency to 1 rad/s for HPF and BANDSTOP filters
29  ******************************************************************************/
30 double DSPL_API filter_ws1(int ord, double rp, double rs, int type)
31 {
32  double es2, ep2, gs2, x, ws;
33 
34  if(ord<1 || rp < 0.0 || rs < 0.0)
35  return -1.0;
36 
37  es2 = pow(10.0, rs*0.1) - 1.0;
38  ep2 = pow(10.0, rp*0.1) - 1.0;
39  gs2 = 1.0 / (1.0 + es2);
40 
41  x = (1.0 - gs2) / (gs2 * ep2);
42 
43  switch( type & DSPL_FILTER_APPROX_MASK)
44  {
45  case DSPL_FILTER_BUTTER:
46  ws = pow(x, 0.5 / (double)ord);
47  break;
48  case DSPL_FILTER_CHEBY1:
49  case DSPL_FILTER_CHEBY2:
50  x = sqrt(x) + sqrt(x - 1.0);
51  x = log(x) / (double)ord;
52  ws = 0.5 * (exp(-x) + exp(x));
53  break;
54  case DSPL_FILTER_ELLIP:
55  {
56  double k, k1;
57  complex_t y, z;
58  int res;
59  k = sqrt(ep2 / es2);
60  res = ellip_modulareq(rp, rs, ord, &k1);
61  if(res != RES_OK)
62  {
63  ws = -1.0;
64  break;
65  }
66  RE(z) = sqrt(x);
67  IM(z) = 0.0;
68 
69  res = ellip_acd_cmplx(&z, 1, k, &y);
70  if(res != RES_OK)
71  {
72  ws = -1.0;
73  break;
74  }
75  RE(y) /= (double)ord;
76  IM(y) /= (double)ord;
77  res = ellip_cd_cmplx(&y, 1, k1, &z);
78  if(res != RES_OK)
79  {
80  ws = -1.0;
81  break;
82  }
83  ws = RE(z);
84  break;
85  }
86  default:
87  ws = -1.0;
88  break;
89  }
90  return ws;
91 }
92 
93 /******************************************************************************
94  * low 2 bandpass transformation
95  ******************************************************************************/
96 int DSPL_API low2bp(double* b, double* a, int ord,
97  double w0, double wpl, double wph,
98  double* beta, double* alpha)
99 {
100 
101  double num[3] = {0.0, 0.0, 1.0};
102  double den[3] = {0.0, 0.0, 0.0};
103 
104  if(!b || !a || !beta || !alpha)
105  return ERROR_PTR;
106  if(ord < 1)
107  return ERROR_FILTER_ORD;
108  if(w0 <= 0.0 || wpl <= 0.0 || wph <= 0.0 || wph <= wpl)
109  return ERROR_FILTER_FT;
110 
111  num[0] = (wph * wpl) / (w0 * w0);
112  den[1] = (wph - wpl) / w0;
113 
114  return ratcompos(b, a, ord, num, den, 2, beta, alpha);
115 }
116 
117 
118 
119 
120 
121 /******************************************************************************
122  * low 2 bandstop transformation
123  ******************************************************************************/
124 int DSPL_API low2bs(double* b, double* a, int ord,
125  double w0, double wsl, double wsh,
126  double* beta, double* alpha)
127 {
128 
129  double den[3] = {0.0, 0.0, 1.0};
130  double num[3] = {0.0, 0.0, 0.0};
131 
132  if(!b || !a || !beta || !alpha)
133  return ERROR_PTR;
134  if(ord < 1)
135  return ERROR_FILTER_ORD;
136  if(w0 <= 0.0 || wsl <= 0.0 || wsh <= 0.0 || wsh <= wsl)
137  return ERROR_FILTER_FT;
138 
139  den[0] = (wsh * wsl) / (w0 * w0);
140  num[1] = (wsh - wsl) / w0;
141 
142  return ratcompos(b, a, ord, num, den, 2, beta, alpha);
143 }
144 
145 
146 
147 
148 /******************************************************************************
149  * low 2 high transformation
150  ******************************************************************************/
151 int DSPL_API low2high(double* b, double* a, int ord, double w0, double w1,
152  double* beta, double* alpha)
153 {
154 
155  double num[2] = {0.0, 0.0};
156  double den[2] = {0.0, 1.0};
157 
158  if(!b || !a || !beta || !alpha)
159  return ERROR_PTR;
160  if(ord < 1)
161  return ERROR_FILTER_ORD;
162  if(w0 <= 0.0 || w1 <= 0.0)
163  return ERROR_FILTER_FT;
164 
165  num[0] = w1 / w0;
166 
167  return ratcompos(b, a, ord, num, den, 1, beta, alpha);
168 }
169 
170 
171 
172 
173 
174 
175 
176 /******************************************************************************
177  low 2 low transformation
178 *******************************************************************************/
179 int DSPL_API low2low(double* b, double* a, int ord, double w0, double w1,
180  double* beta, double* alpha)
181 {
182 
183  double num[2] = {0.0, 1.0};
184  double den[2] = {0.0, 0.0};
185 
186  if(!b || !a || !beta || !alpha)
187  return ERROR_PTR;
188  if(ord < 1)
189  return ERROR_FILTER_ORD;
190  if(w0 <= 0.0 || w1 <= 0.0)
191  return ERROR_FILTER_FT;
192 
193  den[0] = w1 / w0;
194 
195  return ratcompos(b, a, ord, num, den, 1, beta, alpha);
196 }
197 
198 
199 
200 
201 /******************************************************************************
202 Rational composition
203 *******************************************************************************/
204 int DSPL_API ratcompos(double* b, double* a, int n,
205  double* c, double* d, int p,
206  double* beta, double* alpha)
207 {
208 
209  int k2, i, k, pn, pd, ln, ld, k2s, nk2s;
210  double *num = NULL, *den = NULL, *ndn = NULL, *ndd = NULL;
211  int res;
212 
213  if (!a || !b || !c || !d || !beta || !alpha)
214  {
215  res = ERROR_PTR;
216  goto exit_label;
217  }
218  if(n < 1 || p < 1)
219  {
220  res = ERROR_SIZE;
221  goto exit_label;
222  }
223 
224  k2 = (n*p)+1;
225  k2s = k2*sizeof(double); /* alpha and beta size */
226  nk2s = (n+1)*k2*sizeof(double); /* num, den, ndn and ndd size */
227 
228  num = (double*)malloc(nk2s);
229  den = (double*)malloc(nk2s);
230  ndn = (double*)malloc(nk2s);
231  ndd = (double*)malloc(nk2s);
232 
233  memset(num, 0, nk2s);
234  memset(den, 0, nk2s);
235  memset(ndn, 0, nk2s);
236  memset(ndd, 0, nk2s);
237 
238 
239  num[0] = den[0] = 1.0;
240  pn = 0;
241  ln = 1;
242  for(i = 1; i < n+1; i++)
243  {
244  res = conv(num+pn, ln, c, p+1, num+pn+k2);
245  if(res!=RES_OK)
246  goto exit_label;
247  res = conv(den+pn, ln, d, p+1, den+pn+k2);
248  if(res!=RES_OK)
249  goto exit_label;
250  pn += k2;
251  ln += p;
252  }
253 
254  pn = 0;
255  pd = n*k2;
256  ln = 1;
257  ld = k2;
258 
259  for (i = 0; i < n+1; i++)
260  {
261  res = conv(num + pn, ln, den + pd, ld, ndn + i*k2);
262  if(res!=RES_OK)
263  goto exit_label;
264  ln += p;
265  ld -= p;
266  pn += k2;
267  pd -= k2;
268  }
269 
270  for (i = 0; i < n+1; i++)
271  {
272  for (k = 0; k < k2; k++)
273  {
274  ndd[i*k2 + k] = ndn[i*k2 + k] * a[i];
275  ndn[i*k2 + k] *= b[i];
276  }
277  }
278 
279 
280 
281  memset(alpha, 0, k2s);
282  memset(beta, 0, k2s);
283 
284  for (k = 0; k < k2; k++)
285  {
286  for (i = 0; i < n+1; i++)
287  {
288  beta[k] += ndn[i*k2 + k];
289  alpha[k] += ndd[i*k2 + k];
290  }
291  }
292 
293  res = RES_OK;
294 exit_label:
295  if(num)
296  free(num);
297  if(den)
298  free(den);
299  if(ndn)
300  free(ndn);
301  if(ndd)
302  free(ndd);
303 
304  return res;
305 }
306 
int low2low(double *b, double *a, int ord, double w0, double w1, double *beta, double *alpha)
Частотное преобразование ФНЧ-ФНЧ
Definition: filter_ft.c:179
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:148
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:41
#define ERROR_PTR
Ошибка указателя. Данная ошибка означает, что один из обязательных указателей (память под который дол...
Definition: dspl.h:140
int ratcompos(double *b, double *a, int n, double *c, double *d, int p, double *beta, double *alpha)
Рациональная композиця
Definition: filter_ft.c:204
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:78
int low2high(double *b, double *a, int ord, double w0, double w1, double *beta, double *alpha)
Частотное преобразование ФНЧ-ФВЧ
Definition: filter_ft.c:151
#define RE(x)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:77
#define ERROR_FILTER_ORD
Порядок фильтра задан неверно. Порядок фильтра должен быть задан положительным целым значением.
Definition: dspl.h:111
int ellip_cd_cmplx(complex_t *u, int n, double k, complex_t *y)
Эллиптическая функция Якоби комплексного аргумента
Definition: ellipj.c:424
int conv(double *a, int na, double *b, int nb, double *c)
Линейная свертка двух вещественных векторов
Definition: conv.c:84
int ellip_acd_cmplx(complex_t *w, int n, double k, complex_t *u)
Обратная эллиптическая функция Якоби комплексного аргумента
Definition: ellipj.c:125
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:94