mirror
/
libdspl-2.0
kopia lustrzana https://github.com/Dsplib/libdspl-2.0
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

added filter_fir and bessel_i0 function

pull/6/head
Dsplib 2019-02-19 23:49:00 +03:00
rodzic 966d85edb4
commit 242123b597
7 zmienionych plików z 334 dodań i 39 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

3
dspl/src/dspl_internal.h
Wyświetl plik

@ -85,7 +85,8 @@ int win_lanczos (double *w, int n, int win_type);
int win_nuttall (double *w, int n, int win_type);
int win_rect (double *w, int n);
int fir_linphase_lpf(int ord, double wp, int wintype,
double winparam, double* h);
#endif

161
dspl/src/filter_fir.c 100644
Wyświetl plik

@ -0,0 +1,161 @@
/*
* Copyright (c) 2015-2019 Sergey Bakhurin
* Digital Signal Processing Library [http://dsplib.org]
*
* This file is part of libdspl-2.0.
*
* is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* DSPL is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Foobar. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "dspl.h"
#include "dspl_internal.h"
/******************************************************************************
* Linear phase lowpass filter
******************************************************************************/
int fir_linphase_lpf(int ord, double wp, int wintype,
double winparam, double* h)
{
int n, err = RES_OK;
double *w = NULL;
w = (double*)malloc((ord+1)*sizeof(double));
err = linspace(-(double)ord*0.5, (double)ord*0.5, ord+1, DSPL_SYMMETRIC, w);
if(err!=RES_OK)
goto error_proc;
err = sinc(w, ord+1, M_PI*wp, h);
if(err!=RES_OK)
goto error_proc;
err = window(w, ord+1, wintype | DSPL_SYMMETRIC, winparam);
if(err!=RES_OK)
goto error_proc;
for(n = 0; n < ord+1; n++)
h[n] *= w[n] * wp;
error_proc:
if(w)
free(w);
return err;
}
/******************************************************************************
* Linear phase FIR filter
******************************************************************************/
int DSPL_API fir_linphase(int ord, double w0, double w1, int filter_type,
int wintype, double winparam, double* h)
{
int n, err;
double wc, b, del;
if(ord<1)
return ERROR_FILTER_ORD;
if(w0 <= 0.0)
return ERROR_FILTER_WP;
if(!h)
return ERROR_PTR;
switch(filter_type & DSPL_FILTER_TYPE_MASK)
{
// Lowpass FIR coefficients calculation
case DSPL_FILTER_LPF:
err = fir_linphase_lpf(ord, w0, wintype, winparam, h);
break;
// Highpass FIR coefficients calculation
case DSPL_FILTER_HPF:
err = fir_linphase_lpf(ord, 1.0-w0, wintype, winparam, h);
if(err == RES_OK)
{
// LPF filter frequency inversion
for(n = 0; n < ord+1; n+=2)
h[n] = -h[n];
}
break;
// Bandpass FIR coefficients calculation
case DSPL_FILTER_BPASS:
if(w1 < w0)
{
err = ERROR_FILTER_WS;
break;
}
wc = (w0 + w1) * 0.5; // central frequency
b = w1 - w0; // bandwidth
err = fir_linphase_lpf(ord, b*0.5, wintype, winparam, h);
if(err == RES_OK)
{
// LPF frequency shifting to the central frequency
del = 0.5 * (double)ord;
for(n = 0; n < ord+1; n++)
h[n] *= 2.0 * cos(M_PI * ((double)n - del) * wc);
}
break;
// BandStop FIR coefficients calculation
// ATTENTION! Bandstop filter must be even order only!
case DSPL_FILTER_BSTOP:
{
double *h0 = NULL;
// check filter order. Return error if order is odd.
if(ord%2)
return ERROR_FILTER_ORD;
// check frequency (w1 must be higher than w0)
if(w1 < w0)
{
err = ERROR_FILTER_WS;
break;
}
// temp coeff vector
h0 = (double*)malloc((ord+1) * sizeof(double));
// calculate LPF
err = fir_linphase(ord, w0, 0.0, DSPL_FILTER_LPF, wintype, winparam, h0);
if(err!=RES_OK)
{
free(h0);
return err;
}
// calculate HPF
err = fir_linphase(ord, w1, 0.0, DSPL_FILTER_HPF, wintype, winparam, h);
if(err==RES_OK)
{
// Bandstop filter is sum of lowpass nad highpass filters
for(n = 0; n < ord+1; n++)
h[n] += h0[n];
}
free(h0);
break;
}
default:
err = ERROR_FILTER_FT;
}
return err;
}

83
dspl/src/math.c
Wyświetl plik

@ -24,6 +24,89 @@
#include "dspl.h"
/*******************************************************************************
Modified Bessel Function of the First Kind – I0(x) [1]
[1] Rational Approximations for the Modified Bessel Function
of the First Kind – I0(x) for Computations with Double Precision
by PAVEL HOLOBORODKO on NOVEMBER 11, 2015
https://www.advanpix.com/2015/11/11/
*******************************************************************************/
int DSPL_API bessel_i0(double *x, int n, double* y)
{
double P16[17] = { 1.0000000000000000000000801e+00,
2.4999999999999999999629693e-01,
2.7777777777777777805664954e-02,
1.7361111111111110294015271e-03,
6.9444444444444568581891535e-05,
1.9290123456788994104574754e-06,
3.9367598891475388547279760e-08,
6.1511873265092916275099070e-10,
7.5940584360755226536109511e-12,
7.5940582595094190098755663e-14,
6.2760839879536225394314453e-16,
4.3583591008893599099577755e-18,
2.5791926805873898803749321e-20,
1.3141332422663039834197910e-22,
5.9203280572170548134753422e-25,
2.0732014503197852176921968e-27,
1.1497640034400735733456400e-29};
double P22[23] = { 3.9894228040143265335649948e-01,
4.9867785050353992900698488e-02,
2.8050628884163787533196746e-02,
2.9219501690198775910219311e-02,
4.4718622769244715693031735e-02,
9.4085204199017869159183831e-02,
-1.0699095472110916094973951e-01,
2.2725199603010833194037016e+01,
-1.0026890180180668595066918e+03,
3.1275740782277570164423916e+04,
-5.9355022509673600842060002e+05,
2.6092888649549172879282592e+06,
2.3518420447411254516178388e+08,
-8.9270060370015930749184222e+09,
1.8592340458074104721496236e+11,
-2.6632742974569782078420204e+12,
2.7752144774934763122129261e+13,
-2.1323049786724612220362154e+14,
1.1989242681178569338129044e+15,
-4.8049082153027457378879746e+15,
1.3012646806421079076251950e+16,
-2.1363029690365351606041265e+16,
1.6069467093441596329340754e+16};
double x2;
int k;
if(!x || !y)
return ERROR_PTR;
if(n < 1)
return ERROR_SIZE;
for(k =0; k < n; k++)
{
if(x[k] < 0.0)
return ERROR_NEGATIVE;
if(x[k] < 7.75)
{
x2 = x[k] * x[k] * 0.25;
polyval(P16, 16, &x2, 1, y+k);
y[k] = x2 * y[k] + 1.0;
}
else
{
x2 = 1.0 / x[k];
polyval(P22, 22, &x2, 1, y+k);
y[k] *= exp(x[k]) / sqrt(x[k]);
}
}
return RES_OK;
}
/*******************************************************************************
module operator for double
*******************************************************************************/

12
include/dspl.c
Wyświetl plik

@ -38,6 +38,7 @@
p_acos_cmplx acos_cmplx ;
p_asin_cmplx asin_cmplx ;
p_butter_ap butter_ap ;
p_bessel_i0 bessel_i0 ;
p_butter_ap_zp butter_ap_zp ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
@ -74,14 +75,15 @@ p_ifft_cmplx ifft_cmplx ;
p_dspl_info dspl_info ;
p_farrow_lagrange farrow_lagrange ;
p_farrow_spline farrow_spline ;
p_filter_iir filter_iir ;
p_filter_zp2ab filter_zp2ab ;
p_fft fft ;
p_fft_cmplx fft_cmplx ;
p_fft_create fft_create ;
p_fft_free fft_free ;
p_fft_shift fft_shift ;
p_fft_shift_cmplx fft_shift_cmplx ;
p_filter_iir filter_iir ;
p_filter_zp2ab filter_zp2ab ;
p_fir_linphase fir_linphase ;
p_flipip flipip ;
p_flipip_cmplx flipip_cmplx ;
p_fourier_integral_cmplx fourier_integral_cmplx ;
@ -186,6 +188,7 @@ void* dspl_load()
LOAD_FUNC(acos_cmplx);
LOAD_FUNC(asin_cmplx);
LOAD_FUNC(bessel_i0);
LOAD_FUNC(butter_ap);
LOAD_FUNC(butter_ap_zp);
LOAD_FUNC(cheby_poly1);
@ -223,14 +226,15 @@ void* dspl_load()
LOAD_FUNC(dspl_info);
LOAD_FUNC(farrow_lagrange);
LOAD_FUNC(farrow_spline);
LOAD_FUNC(filter_iir);
LOAD_FUNC(filter_zp2ab);
LOAD_FUNC(fft);
LOAD_FUNC(fft_cmplx);
LOAD_FUNC(fft_create);
LOAD_FUNC(fft_free);
LOAD_FUNC(fft_shift);
LOAD_FUNC(fft_shift_cmplx);
LOAD_FUNC(filter_iir);
LOAD_FUNC(filter_zp2ab);
LOAD_FUNC(fir_linphase);
LOAD_FUNC(flipip);
LOAD_FUNC(flipip_cmplx);
LOAD_FUNC(fourier_integral_cmplx);

51
include/dspl.h
Wyświetl plik

@ -97,6 +97,8 @@ typedef struct
#define ERROR_FILTER_ORD 0x06101518
#define ERROR_FILTER_RP 0x06101816
#define ERROR_FILTER_RS 0x06101819
#define ERROR_FILTER_WP 0x06102316
#define ERROR_FILTER_WS 0x06102319
#define ERROR_FNAME 0x06140113
#define ERROR_FOPEN 0x06151605
#define ERROR_FREAD_SIZE 0x06180501
@ -172,6 +174,13 @@ typedef struct
#define DSPL_WIN_COS 0x00040000
#define DSPL_WIN_CHEBY 0x00080000
#define DSPL_FILTER_TYPE_MASK 0x0000000F
#define DSPL_FILTER_LPF 0x00000001
#define DSPL_FILTER_HPF 0x00000002
#define DSPL_FILTER_BPASS 0x00000004
#define DSPL_FILTER_BSTOP 0x00000008
#define ELLIP_ITER 16
#define ELLIP_MAX_ORD 24
@ -221,6 +230,10 @@ DECLARE_FUNC(int, asin_cmplx, complex_t*
COMMA int
COMMA complex_t*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, bessel_i0, double* x
COMMA int n
COMMA double* y);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, butter_ap, double
COMMA int
COMMA double*
@ -420,21 +433,6 @@ DECLARE_FUNC(int, farrow_spline, double*
COMMA double**
COMMA int*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, filter_iir, double*
COMMA double*
COMMA int
COMMA double*
COMMA int
COMMA double*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, filter_zp2ab, complex_t*
COMMA int
COMMA complex_t*
COMMA int
COMMA int
COMMA double*
COMMA double*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fft, double*
COMMA int
COMMA fft_t*
@ -458,6 +456,29 @@ DECLARE_FUNC(int, fft_shift_cmplx, complex_t*
COMMA int
COMMA complex_t*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, filter_iir, double*
COMMA double*
COMMA int
COMMA double*
COMMA int
COMMA double*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, filter_zp2ab, complex_t*
COMMA int
COMMA complex_t*
COMMA int
COMMA int
COMMA double*
COMMA double*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fir_linphase, int ord
COMMA double w0
COMMA double w1
COMMA int filter_type
COMMA int wintype
COMMA double winparam
COMMA double* h);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, flipip, double*
COMMA int);
//------------------------------------------------------------------------------

12
release/include/dspl.c
Wyświetl plik

@ -38,6 +38,7 @@
p_acos_cmplx acos_cmplx ;
p_asin_cmplx asin_cmplx ;
p_butter_ap butter_ap ;
p_bessel_i0 bessel_i0 ;
p_butter_ap_zp butter_ap_zp ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
@ -74,14 +75,15 @@ p_ifft_cmplx ifft_cmplx ;
p_dspl_info dspl_info ;
p_farrow_lagrange farrow_lagrange ;
p_farrow_spline farrow_spline ;
p_filter_iir filter_iir ;
p_filter_zp2ab filter_zp2ab ;
p_fft fft ;
p_fft_cmplx fft_cmplx ;
p_fft_create fft_create ;
p_fft_free fft_free ;
p_fft_shift fft_shift ;
p_fft_shift_cmplx fft_shift_cmplx ;
p_filter_iir filter_iir ;
p_filter_zp2ab filter_zp2ab ;
p_fir_linphase fir_linphase ;
p_flipip flipip ;
p_flipip_cmplx flipip_cmplx ;
p_fourier_integral_cmplx fourier_integral_cmplx ;
@ -186,6 +188,7 @@ void* dspl_load()
LOAD_FUNC(acos_cmplx);
LOAD_FUNC(asin_cmplx);
LOAD_FUNC(bessel_i0);
LOAD_FUNC(butter_ap);
LOAD_FUNC(butter_ap_zp);
LOAD_FUNC(cheby_poly1);
@ -223,14 +226,15 @@ void* dspl_load()
LOAD_FUNC(dspl_info);
LOAD_FUNC(farrow_lagrange);
LOAD_FUNC(farrow_spline);
LOAD_FUNC(filter_iir);
LOAD_FUNC(filter_zp2ab);
LOAD_FUNC(fft);
LOAD_FUNC(fft_cmplx);
LOAD_FUNC(fft_create);
LOAD_FUNC(fft_free);
LOAD_FUNC(fft_shift);
LOAD_FUNC(fft_shift_cmplx);
LOAD_FUNC(filter_iir);
LOAD_FUNC(filter_zp2ab);
LOAD_FUNC(fir_linphase);
LOAD_FUNC(flipip);
LOAD_FUNC(flipip_cmplx);
LOAD_FUNC(fourier_integral_cmplx);

51
release/include/dspl.h
Wyświetl plik

@ -97,6 +97,8 @@ typedef struct
#define ERROR_FILTER_ORD 0x06101518
#define ERROR_FILTER_RP 0x06101816
#define ERROR_FILTER_RS 0x06101819
#define ERROR_FILTER_WP 0x06102316
#define ERROR_FILTER_WS 0x06102319
#define ERROR_FNAME 0x06140113
#define ERROR_FOPEN 0x06151605
#define ERROR_FREAD_SIZE 0x06180501
@ -172,6 +174,13 @@ typedef struct
#define DSPL_WIN_COS 0x00040000
#define DSPL_WIN_CHEBY 0x00080000
#define DSPL_FILTER_TYPE_MASK 0x0000000F
#define DSPL_FILTER_LPF 0x00000001
#define DSPL_FILTER_HPF 0x00000002
#define DSPL_FILTER_BPASS 0x00000004
#define DSPL_FILTER_BSTOP 0x00000008
#define ELLIP_ITER 16
#define ELLIP_MAX_ORD 24
@ -221,6 +230,10 @@ DECLARE_FUNC(int, asin_cmplx, complex_t*
COMMA int
COMMA complex_t*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, bessel_i0, double* x
COMMA int n
COMMA double* y);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, butter_ap, double
COMMA int
COMMA double*
@ -420,21 +433,6 @@ DECLARE_FUNC(int, farrow_spline, double*
COMMA double**
COMMA int*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, filter_iir, double*
COMMA double*
COMMA int
COMMA double*
COMMA int
COMMA double*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, filter_zp2ab, complex_t*
COMMA int
COMMA complex_t*
COMMA int
COMMA int
COMMA double*
COMMA double*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fft, double*
COMMA int
COMMA fft_t*
@ -458,6 +456,29 @@ DECLARE_FUNC(int, fft_shift_cmplx, complex_t*
COMMA int
COMMA complex_t*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, filter_iir, double*
COMMA double*
COMMA int
COMMA double*
COMMA int
COMMA double*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, filter_zp2ab, complex_t*
COMMA int
COMMA complex_t*
COMMA int
COMMA int
COMMA double*
COMMA double*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fir_linphase, int ord
COMMA double w0
COMMA double w1
COMMA int filter_type
COMMA int wintype
COMMA double winparam
COMMA double* h);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, flipip, double*
COMMA int);
//------------------------------------------------------------------------------