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 cheby1 filter

pull/2/head
Dsplib 2018-04-03 23:15:14 +03:00
rodzic 5dfdb06091
commit 3d029712f3
12 zmienionych plików z 191 dodań i 14 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

137
dspl/src/filter_ap.c
Wyświetl plik

@ -30,14 +30,14 @@
/**************************************************************************************************
Analog Normalized Butterworth filter
***************************************************************************************************/
int DSPL_API butter_ap(double Rp, int ord, double* b, double* a)
int DSPL_API butter_ap(double rp, int ord, double* b, double* a)
{
int res;
int res;
complex_t *z = NULL;
complex_t *p = NULL;
int nz, np;
if(Rp < 0.0)
if(rp < 0.0)
return ERROR_FILTER_RP;
if(ord < 1)
return ERROR_FILTER_ORD;
@ -48,7 +48,7 @@ int DSPL_API butter_ap(double Rp, int ord, double* b, double* a)
p = (complex_t*) malloc(ord*sizeof(complex_t));
res = butter_ap_zp(ord, Rp, z, &nz, p, &np);
res = butter_ap_zp(ord, rp, z, &nz, p, &np);
if(res != RES_OK)
goto exit_label;
@ -76,23 +76,23 @@ Analog Normalized Butterworth filter zeros and poles
***************************************************************************************************/
int DSPL_API butter_ap_zp(int ord, double rp, complex_t *z, int* nz, complex_t *p, int* np)
{
double alpha;
double theta;
double alpha;
double theta;
double ep;
int r;
int L;
int L;
int ind = 0, k;
if(rp < 0 || rp == 0)
return ERROR_FILTER_RP;
if(ord < 1)
return ERROR_FILTER_ORD;
if(!z || !p || !nz || !np)
return ERROR_PTR;
return ERROR_FILTER_RP;
if(ord < 1)
return ERROR_FILTER_ORD;
if(!z || !p || !nz || !np)
return ERROR_PTR;
ep = sqrt(pow(10.0, rp*0.1) - 1.0);
r = ord % 2;
L = (int)((ord-r)/2);
r = ord % 2;
L = (int)((ord-r)/2);
alpha = pow(ep, -1.0/(double)ord);
if(r)
@ -119,6 +119,115 @@ int DSPL_API butter_ap_zp(int ord, double rp, complex_t *z, int* nz, complex_t *
/**************************************************************************************************
Analog Normalized Chebyshev type 1 filter
***************************************************************************************************/
int DSPL_API cheby1_ap(double rp, int ord, double* b, double* a)
{
int res;
complex_t *z = NULL;
complex_t *p = NULL;
int nz, np, k;
complex_t h0 = {1.0, 0.0};
double tmp;
if(rp < 0.0)
return ERROR_FILTER_RP;
if(ord < 1)
return ERROR_FILTER_ORD;
if(!a || !b)
return ERROR_PTR;
z = (complex_t*) malloc(ord*sizeof(complex_t));
p = (complex_t*) malloc(ord*sizeof(complex_t));
res = cheby1_ap_zp(ord, rp, z, &nz, p, &np);
if(res != RES_OK)
goto exit_label;
res = filter_zp2ab(z, nz, p, np, ord, b, a);
if(res != RES_OK)
goto exit_label;
if(!(ord % 2))
RE(h0) = 1.0 / pow(10.0, rp*0.05);
for(k = 0; k < np; k++)
{
tmp = CMRE(h0, p[k]);
IM(h0) = CMIM(h0, p[k]);
RE(h0) = tmp;
}
b[0] = fabs(RE(h0));
exit_label:
if(z)
free(z);
if(p)
free(p);
return res;
}
/**************************************************************************************************
Analog Normalized Chebyshev type 1 filter zeros and poles
***************************************************************************************************/
int DSPL_API cheby1_ap_zp(int ord, double rp, complex_t *z, int* nz, complex_t *p, int* np)
{
double alpha;
double theta;
double ep;
double beta;
double shbeta;
double chbeta;
int r;
int L;
int ind = 0, k;
if(rp < 0 || rp == 0)
return ERROR_FILTER_RP;
if(ord < 1)
return ERROR_FILTER_ORD;
if(!z || !p || !nz || !np)
return ERROR_PTR;
ep = sqrt(pow(10.0, rp*0.1) - 1.0);
r = ord % 2;
L = (int)((ord-r)/2);
beta = asinh(1.0/ep)/(double)ord;
chbeta = cosh(beta);
shbeta = sinh(beta);
if(r)
{
RE(p[ind]) = -shbeta;
IM(p[ind]) = 0.0;
ind++;
}
for(k = 0; k < L; k++)
{
theta = M_PI*(double)(2*k + 1)/(double)(2*ord);
RE(p[ind]) = RE(p[ind+1]) = -shbeta * sin(theta);
IM(p[ind]) = chbeta * cos(theta);
IM(p[ind+1]) = -IM(p[ind]);
ind+=2;
}
*np = ord;
*nz = 0;
return RES_OK;
}
/**************************************************************************************************

4
include/dspl.c
Wyświetl plik

@ -41,6 +41,8 @@ p_butter_ap butter_ap ;
p_butter_ap_zp butter_ap_zp ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
p_cheby1_ap cheby1_ap ;
p_cheby1_ap_zp cheby1_ap_zp ;
p_cmplx2re cmplx2re ;
p_concat concat ;
p_conv conv ;
@ -139,6 +141,8 @@ void* dspl_load()
LOAD_FUNC(butter_ap_zp);
LOAD_FUNC(cheby_poly1);
LOAD_FUNC(cheby_poly2);
LOAD_FUNC(cheby1_ap);
LOAD_FUNC(cheby1_ap_zp);
LOAD_FUNC(cmplx2re);
LOAD_FUNC(concat);
LOAD_FUNC(conv);

2
include/dspl.h
Wyświetl plik

@ -200,6 +200,8 @@ DECLARE_FUNC(int, butter_ap_zp, int COMMA double COMMA comple
DECLARE_FUNC(int, cheby_poly1, double* COMMA int COMMA int COMMA double*);
DECLARE_FUNC(int, cheby_poly2, double* COMMA int COMMA int COMMA double*);
DECLARE_FUNC(int, cheby1_ap, double COMMA int COMMA double* COMMA double*);
DECLARE_FUNC(int, cheby1_ap_zp, int COMMA double COMMA complex_t* COMMA int* COMMA complex_t* COMMA int*);
DECLARE_FUNC(int, cmplx2re, complex_t* COMMA int COMMA double* COMMA double*);
DECLARE_FUNC(int, concat, void* COMMA size_t COMMA void* COMMA size_t COMMA void*);
DECLARE_FUNC(int, conv, double* COMMA int COMMA double* COMMA int COMMA double*);

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

@ -41,6 +41,8 @@ p_butter_ap butter_ap ;
p_butter_ap_zp butter_ap_zp ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
p_cheby1_ap cheby1_ap ;
p_cheby1_ap_zp cheby1_ap_zp ;
p_cmplx2re cmplx2re ;
p_concat concat ;
p_conv conv ;
@ -139,6 +141,8 @@ void* dspl_load()
LOAD_FUNC(butter_ap_zp);
LOAD_FUNC(cheby_poly1);
LOAD_FUNC(cheby_poly2);
LOAD_FUNC(cheby1_ap);
LOAD_FUNC(cheby1_ap_zp);
LOAD_FUNC(cmplx2re);
LOAD_FUNC(concat);
LOAD_FUNC(conv);

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

@ -200,6 +200,8 @@ DECLARE_FUNC(int, butter_ap_zp, int COMMA double COMMA comple
DECLARE_FUNC(int, cheby_poly1, double* COMMA int COMMA int COMMA double*);
DECLARE_FUNC(int, cheby_poly2, double* COMMA int COMMA int COMMA double*);
DECLARE_FUNC(int, cheby1_ap, double COMMA int COMMA double* COMMA double*);
DECLARE_FUNC(int, cheby1_ap_zp, int COMMA double COMMA complex_t* COMMA int* COMMA complex_t* COMMA int*);
DECLARE_FUNC(int, cmplx2re, complex_t* COMMA int COMMA double* COMMA double*);
DECLARE_FUNC(int, concat, void* COMMA size_t COMMA void* COMMA size_t COMMA void*);
DECLARE_FUNC(int, conv, double* COMMA int COMMA double* COMMA int COMMA double*);

3
test/bin/.gitignore vendored
Wyświetl plik

@ -0,0 +1,3 @@
*.exe
*.so
*.dll

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

Plik binarny nie jest wyświetlany.

1
test/bin/gnuplot/butter_ap_test.plt
Wyświetl plik

@ -4,4 +4,5 @@ unset key
set grid
set xlabel "frequency, rad/s"
set ylabel "Butterworth filter magnitude, dB"
set yrange [-100:5]
plot 'dat/butter_ap_test_mag.txt' with lines

8
test/bin/gnuplot/cheby1_ap_test.plt 100644
Wyświetl plik

@ -0,0 +1,8 @@
set logscale x
unset key
set grid
set xlabel "frequency, rad/s"
set ylabel "Chebyshev type 1 filter magnitude, dB"
set yrange [-100:5]
plot 'dat/cheby1_ap_test_mag.txt' with lines

BIN
test/bin/img/butter_ap_test_mag.png
Wyświetl plik

Plik binarny nie jest wyświetlany.
Side by Side Swipe Nakładka

Przed

Szerokość:  |  Wysokość:  |  Rozmiar: 23 KiB

Po

Szerokość:  |  Wysokość:  |  Rozmiar: 21 KiB

BIN
test/bin/img/cheby1_ap_test_mag.png 100644
Wyświetl plik

Plik binarny nie jest wyświetlany.
Side by Side

Po

Szerokość:  |  Wysokość:  |  Rozmiar: 22 KiB

44
test/src/cheby1_ap_test.c 100644
Wyświetl plik

@ -0,0 +1,44 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define ORD 4
#define N 1000
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
double a[ORD+1], b[ORD+1];
double Rp = 3.0;
int k;
int res = cheby1_ap(Rp, ORD, b, a);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
for(k = 0; k < ORD+1; k++)
printf("b[%2d] = %9.3f a[%2d] = %9.3f\n", k, b[k], k, a[k]);
double w[N], hdb[N];
complex_t h[N];
logspace(-2.0, 2.0, N , DSPL_SYMMETRIC, w);
freqs(b, a, ORD, w, N, h);
for(k = 0; k < N; k++)
hdb[k] = 10.0 * log10(ABSSQR(h[k]));
writetxt(w, hdb, N, "dat/cheby1_ap_test_mag.txt");
dspl_free(handle); // free dspl handle
system("gnuplot -p gnuplot/cheby1_ap_test.plt");
return 0;
}