libdspl-2.0
Digital Signal Processing Algorithm Library
filter_fir.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  * Linear phase lowpass filter
30  ******************************************************************************/
31 int fir_linphase_lpf(int ord, double wp, int win_type,
32  double win_param, double* h)
33 {
34  int n, err = RES_OK;
35  double *w = NULL;
36 
37 
38  w = (double*)malloc((ord+1)*sizeof(double));
39 
40  err = linspace(-(double)ord*0.5, (double)ord*0.5, ord+1, DSPL_SYMMETRIC, w);
41 
42  if(err!=RES_OK)
43  goto error_proc;
44 
45  err = sinc(w, ord+1, M_PI*wp, h);
46 
47  if(err!=RES_OK)
48  goto error_proc;
49 
50  err = window(w, ord+1, win_type | DSPL_SYMMETRIC, win_param);
51 
52  if(err!=RES_OK)
53  goto error_proc;
54 
55  for(n = 0; n < ord+1; n++)
56  h[n] *= w[n] * wp;
57 
58 error_proc:
59  if(w)
60  free(w);
61  return err;
62 }
63 
64 
65 
66 
67 
68 #ifdef DOXYGEN_ENGLISH
69 
190 #endif
191 #ifdef DOXYGEN_RUSSIAN
192 
320 #endif
321 int DSPL_API fir_linphase(int ord, double w0, double w1, int filter_type,
322  int win_type, double win_param, double* h)
323 {
324  int n, err;
325  double wc, b, del;
326 
327  if(ord<1)
328  return ERROR_FILTER_ORD;
329  if(w0 <= 0.0)
330  return ERROR_FILTER_WP;
331  if(!h)
332  return ERROR_PTR;
333 
334  switch(filter_type & DSPL_FILTER_TYPE_MASK)
335  {
336  /* Lowpass FIR coefficients calculation */
337  case DSPL_FILTER_LPF:
338  err = fir_linphase_lpf(ord, w0, win_type, win_param, h);
339  break;
340 
341  /* Highpass FIR coefficients calculation */
342  case DSPL_FILTER_HPF:
343  err = fir_linphase_lpf(ord, 1.0-w0, win_type, win_param, h);
344  if(err == RES_OK)
345  {
346  /* LPF filter frequency inversion */
347  for(n = 0; n < ord+1; n+=2)
348  h[n] = -h[n];
349  }
350  break;
351 
352  /* Bandpass FIR coefficients calculation */
353  case DSPL_FILTER_BPASS:
354  if(w1 < w0)
355  {
356  err = ERROR_FILTER_WS;
357  break;
358  }
359  wc = (w0 + w1) * 0.5; /* central frequency */
360  b = w1 - w0; /* bandwidth */
361  err = fir_linphase_lpf(ord, b*0.5, win_type, win_param, h);
362  if(err == RES_OK)
363  {
364  /* LPF frequency shifting to the central frequency */
365  del = 0.5 * (double)ord;
366  for(n = 0; n < ord+1; n++)
367  h[n] *= 2.0 * cos(M_PI * ((double)n - del) * wc);
368  }
369  break;
370 
371  /* BandStop FIR coefficients calculation */
372  /* ATTENTION! Bandstop filter must be even order only! */
373  case DSPL_FILTER_BSTOP:
374  {
375  double *h0 = NULL;
376 
377  /* check filter order. Return error if order is odd. */
378  if(ord%2)
379  return ERROR_FILTER_ORD;
380 
381  /* check frequency (w1 must be higher than w0) */
382  if(w1 < w0)
383  {
384  err = ERROR_FILTER_WS;
385  break;
386  }
387  /* temp coeff vector */
388  h0 = (double*)malloc((ord+1) * sizeof(double));
389 
390  /* calculate LPF */
391  err = fir_linphase(ord, w0, 0.0, DSPL_FILTER_LPF,
392  win_type, win_param, h0);
393  if(err!=RES_OK)
394  {
395  free(h0);
396  return err;
397  }
398  /* calculate HPF */
399  err = fir_linphase(ord, w1, 0.0, DSPL_FILTER_HPF,
400  win_type, win_param, h);
401  if(err==RES_OK)
402  {
403  /* Bandstop filter is sum of lowpass and highpass filters */
404  for(n = 0; n < ord+1; n++)
405  h[n] += h0[n];
406  }
407  free(h0);
408  break;
409  }
410  default:
411  err = ERROR_FILTER_FT;
412  }
413  return err;
414 }
int DSPL_API window(double *w, int n, int win_type, double param)
Window function calculation.
Definition: win.c:308
#define ERROR_FILTER_WS
The filter frequency parameter is set incorrectly. The frequency must be a positive from 0 to 1.
Definition: dspl.h:520
#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 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_FILTER_FT
The conversion frequencies of the low-pass filter and low-pass filter are incorrect....
Definition: dspl.h:513
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:497
#define ERROR_FILTER_WP
The filter cutoff frequency parameter is set incorrectly. The cutoff frequency of the filter must be ...
Definition: dspl.h:519
int DSPL_API linspace(double x0, double x1, int n, int type, double *x)
Function fills a vector with n linearly spaced elements between x0 and x1.
Definition: array.c:1009
int DSPL_API sinc(double *x, int n, double a, double *y)
Function for the real vector x.
Definition: math.c:936
int DSPL_API fir_linphase(int ord, double w0, double w1, int filter_type, int win_type, double win_param, double *h)
Function calculates linear-phase FIR filter coefficients by window method.
Definition: filter_fir.c:321