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 freqz freqs

pull/2/head
Dsplib 2018-04-01 23:08:42 +03:00
rodzic 6fb251d989
commit 032d4cba74
6 zmienionych plików z 191 dodań i 8 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

171
dspl/src/filter_an.c 100644
Wyświetl plik

@ -0,0 +1,171 @@
/*
* Copyright (c) 2015-2018 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 <string.h>
#include <math.h>
#include "dspl.h"
/**************************************************************************************************
Complex frequency response of an analog filter H(s)
***************************************************************************************************/
int DSPL_API freqs(double* b, double* a, int ord, double* w, int n, complex_t *h)
{
complex_t jw;
complex_t *bc = NULL;
complex_t *ac = NULL;
complex_t num, den;
double mag;
int k;
int res;
if(!b || !a || !w || !h)
return ERROR_PTR;
if(ord<0)
return ERROR_FILTER_ORD;
if(n<1)
return ERROR_SIZE;
RE(jw) = 0.0;
bc = (complex_t*) malloc((ord+1) * sizeof(complex_t));
res = re2cmplx(b, ord+1, bc);
if( res!=RES_OK )
goto exit_label;
ac = (complex_t*) malloc((ord+1) * sizeof(complex_t));
res = re2cmplx(a, ord+1, ac);
if( res!=RES_OK )
goto exit_label;
for(k = 0; k < n; k++)
{
IM(jw) = w[k];
res = polyval_cmplx(bc, ord, &jw, 1, &num);
if(res != RES_OK)
goto exit_label;
res = polyval_cmplx(ac, ord, &jw, 1, &den);
if(res != RES_OK)
goto exit_label;
mag = ABSSQR(den);
if(mag == 0.0)
{
res = ERROR_DIV_ZERO;
goto exit_label;
}
mag = 1.0 / mag;
RE(h[k]) = CMCONJRE(num, den) * mag;
IM(h[k]) = CMCONJIM(num, den) * mag;
}
res = RES_OK;
exit_label:
if(bc)
free(bc);
if(ac)
free(ac);
return res;
}
/**************************************************************************************************
Complex frequency response of a digital filter H(z)
**************************************************************************************************/
int DSPL_API freqz(double* b, double* a, int ord, double* w, int n, complex_t *h)
{
complex_t jw;
complex_t *bc = NULL;
complex_t *ac = NULL;
complex_t num, den;
double mag;
int k;
int res;
if(!b || !w || !h)
return ERROR_PTR;
if(ord<0)
return ERROR_FILTER_ORD;
if(n<1)
return ERROR_SIZE;
bc = (complex_t*) malloc((ord+1) * sizeof(complex_t));
res = re2cmplx(b, ord+1, bc);
if( res!=RES_OK )
goto exit_label;
if(a)
{
// IIR filter if a != NULL
ac = (complex_t*) malloc((ord+1) * sizeof(complex_t));
res = re2cmplx(a, ord+1, ac);
if( res!=RES_OK )
goto exit_label;
for(k = 0; k < n; k++)
{
RE(jw) = cos(w[k]);
IM(jw) = -sin(w[k]);
res = polyval_cmplx(bc, ord, &jw, 1, &num);
if(res != RES_OK)
goto exit_label;
res = polyval_cmplx(ac, ord, &jw, 1, &den);
if(res != RES_OK)
goto exit_label;
mag = ABSSQR(den);
if(mag == 0.0)
{
res = ERROR_DIV_ZERO;
goto exit_label;
}
mag = 1.0 / mag;
RE(h[k]) = CMCONJRE(num, den) * mag;
IM(h[k]) = CMCONJIM(num, den) * mag;
}
}
else
{
// FIR filter if a == NULL
for(k = 0; k < n; k++)
{
RE(jw) = cos(w[k]);
IM(jw) = -sin(w[k]);
res = polyval_cmplx(bc, ord, &jw, 1, h+k);
if(res != RES_OK)
goto exit_label;
}
}
res = RES_OK;
exit_label:
if(bc)
free(bc);
if(ac)
free(ac);
return res;
}

12
include/dspl.c
Wyświetl plik

@ -39,15 +39,15 @@
p_acos_cmplx acos_cmplx ;
p_asin_cmplx asin_cmplx ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
p_cmplx2re cmplx2re ;
p_concat concat ;
p_conv conv ;
p_conv_cmplx conv_cmplx ;
p_cos_cmplx cos_cmplx ;
p_dft dft ;
p_dft_cmplx dft_cmplx ;
p_dft dft ;
p_dft_cmplx dft_cmplx ;
p_dmod dmod ;
p_farrow_lagrange farrow_lagrange ;
p_farrow_spline farrow_spline ;
@ -62,6 +62,8 @@ p_flipip flipip ;
p_flipip_cmplx flipip_cmplx ;
p_fourier_series_dec fourier_series_dec ;
p_fourier_series_rec fourier_series_rec ;
p_freqs freqs ;
p_freqz freqz ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
p_linspace linspace ;
@ -156,6 +158,8 @@ void* dspl_load()
LOAD_FUNC(flipip_cmplx);
LOAD_FUNC(fourier_series_dec);
LOAD_FUNC(fourier_series_rec);
LOAD_FUNC(freqz);
LOAD_FUNC(freqs);
LOAD_FUNC(goertzel);
LOAD_FUNC(goertzel_cmplx);
LOAD_FUNC(linspace);

2
include/dspl.h
Wyświetl plik

@ -218,6 +218,8 @@ DECLARE_FUNC(int, flipip, double* COMMA int);
DECLARE_FUNC(int, flipip_cmplx, complex_t* COMMA int);
DECLARE_FUNC(int, fourier_series_dec, double* COMMA double* COMMA int COMMA double COMMA int COMMA double* COMMA complex_t*);
DECLARE_FUNC(int, fourier_series_rec, double* COMMA complex_t* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, freqs, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, freqz, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel, double* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel_cmplx, complex_t* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, linspace, double COMMA double COMMA int COMMA int COMMA double*);

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

@ -39,15 +39,15 @@
p_acos_cmplx acos_cmplx ;
p_asin_cmplx asin_cmplx ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
p_cmplx2re cmplx2re ;
p_concat concat ;
p_conv conv ;
p_conv_cmplx conv_cmplx ;
p_cos_cmplx cos_cmplx ;
p_dft dft ;
p_dft_cmplx dft_cmplx ;
p_dft dft ;
p_dft_cmplx dft_cmplx ;
p_dmod dmod ;
p_farrow_lagrange farrow_lagrange ;
p_farrow_spline farrow_spline ;
@ -62,6 +62,8 @@ p_flipip flipip ;
p_flipip_cmplx flipip_cmplx ;
p_fourier_series_dec fourier_series_dec ;
p_fourier_series_rec fourier_series_rec ;
p_freqs freqs ;
p_freqz freqz ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
p_linspace linspace ;
@ -156,6 +158,8 @@ void* dspl_load()
LOAD_FUNC(flipip_cmplx);
LOAD_FUNC(fourier_series_dec);
LOAD_FUNC(fourier_series_rec);
LOAD_FUNC(freqz);
LOAD_FUNC(freqs);
LOAD_FUNC(goertzel);
LOAD_FUNC(goertzel_cmplx);
LOAD_FUNC(linspace);

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

@ -218,6 +218,8 @@ DECLARE_FUNC(int, flipip, double* COMMA int);
DECLARE_FUNC(int, flipip_cmplx, complex_t* COMMA int);
DECLARE_FUNC(int, fourier_series_dec, double* COMMA double* COMMA int COMMA double COMMA int COMMA double* COMMA complex_t*);
DECLARE_FUNC(int, fourier_series_rec, double* COMMA complex_t* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, freqs, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, freqz, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel, double* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel_cmplx, complex_t* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, linspace, double COMMA double COMMA int COMMA int COMMA double*);

BIN
test/bin/dft_test.exe
Wyświetl plik

Plik binarny nie jest wyświetlany.