From 4306ab98ee2b3fd11a70f495eaff7592f534f241 Mon Sep 17 00:00:00 2001 From: Dsplib Date: Wed, 23 May 2018 23:36:00 +0300 Subject: [PATCH] fixed bugs in ellip_asn_cmplx and ellip_acd_cmplx. Added ellip_ap_zp --- dspl/src/ellipj.c | 8 ++-- dspl/src/filter_ap.c | 102 +++++++++++++++++++++++++++++++++++++++++ include/dspl.c | 2 + include/dspl.h | 8 ++++ release/include/dspl.c | 2 + release/include/dspl.h | 8 ++++ 6 files changed, 126 insertions(+), 4 deletions(-) diff --git a/dspl/src/ellipj.c b/dspl/src/ellipj.c index ed7d725..65f653a 100644 --- a/dspl/src/ellipj.c +++ b/dspl/src/ellipj.c @@ -101,8 +101,8 @@ int DSPL_API ellip_acd_cmplx(complex_t* w, int n, double k, complex_t* u) t = 2.0 / ABSSQR(tmp1); - RE(tmp0) = t * CMCONJRE(tmp1, u[m]); - IM(tmp0) = t * CMCONJIM(tmp1, u[m]); + RE(tmp0) = t * CMCONJRE(u[m], tmp1); + IM(tmp0) = t * CMCONJIM(u[m], tmp1); RE(u[m]) = RE(tmp0); IM(u[m]) = IM(tmp0); @@ -195,8 +195,8 @@ int DSPL_API ellip_asn_cmplx(complex_t* w, int n, double k, complex_t* u) t = 2.0 / ABSSQR(tmp1); - RE(tmp0) = t * CMCONJRE(tmp1, u[m]); - IM(tmp0) = t * CMCONJIM(tmp1, u[m]); + RE(tmp0) = t * CMCONJRE(u[m], tmp1); + IM(tmp0) = t * CMCONJIM(u[m], tmp1); RE(u[m]) = RE(tmp0); IM(u[m]) = IM(tmp0); diff --git a/dspl/src/filter_ap.c b/dspl/src/filter_ap.c index 1f6aca2..ca8762c 100644 --- a/dspl/src/filter_ap.c +++ b/dspl/src/filter_ap.c @@ -351,6 +351,108 @@ int DSPL_API cheby2_ap_zp(int ord, double rs, complex_t *z, int* nz, + + + + + +/****************************************************************************** + A *nalog Normalized Chebyshev type 2 filter zeros and poles + *******************************************************************************/ +int DSPL_API ellip_ap_zp(int ord, double rp, double rs, + complex_t *z, int* nz, complex_t *p, int* np) +{ + double es, ep; + int L, r, n, res; + int iz, ip; + double ke, k, u, t; + complex_t tc, v0, jv0; + + + if(rp < 0 || rp == 0) + return ERROR_FILTER_RP; + if(rs < 0 || rs == 0) + return ERROR_FILTER_RS; + if(ord < 1) + return ERROR_FILTER_ORD; + if(!z || !p || !nz || !np) + return ERROR_PTR; + + es = sqrt(pow(10.0, rs*0.1) - 1.0); + ep = sqrt(pow(10.0, rp*0.1) - 1.0); + ke = ep / es; + + r = ord % 2; + L = (int)((ord-r)/2); + + res = ellip_modulareq(rp, rs, ord, &k); + if(res != RES_OK) + return res; + // v0 + RE(tc) = 0.0; + IM(tc) = 1.0 / ep; + + ellip_asn_cmplx(&tc, 1, ke, &v0); + + t = RE(v0); + RE(v0) = IM(v0) / (double)ord; + IM(v0) = -t / (double)ord; + + RE(jv0) = -IM(v0); + IM(jv0) = RE(v0); + + + iz = ip = 0; + + if(r) + { + res = ellip_sn_cmplx(&jv0, 1, k, &tc); + if(res != RES_OK) + return res; + RE(p[0]) = -IM(tc); + IM(p[0]) = RE(tc); + ip = 1; + } + + for(n = 0; n < L; n++) + { + u = (double)(2 * n + 1)/(double)ord; + + res = ellip_cd(& u, 1, k, &t); + if(res != RES_OK) + return res; + + RE(z[iz]) = RE(z[iz+1]) = 0.0; + IM(z[iz]) = 1.0/(k*t); + IM(z[iz+1]) = -1.0/(k*t); + iz+=2; + + RE(tc) = u - RE(jv0); + IM(tc) = - IM(jv0); + + res = ellip_cd_cmplx(&tc, 1, k, p+ip+1); + if(res != RES_OK) + return res; + + RE(p[ip]) = -IM(p[ip+1]); + IM(p[ip]) = RE(p[ip+1]); + + RE(p[ip+1]) = RE(p[ip]); + IM(p[ip+1]) = -IM(p[ip]); + + ip+=2; + + + } + *nz = iz; + *np = ip; + + return RES_OK; +} + + + + /****************************************************************************** Zeros and poles to filter coefficients recalc *******************************************************************************/ diff --git a/include/dspl.c b/include/dspl.c index a3fc94e..c4e57e4 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -56,6 +56,7 @@ p_dft_cmplx dft_cmplx ; p_dmod dmod ; p_ellip_acd ellip_acd ; p_ellip_acd_cmplx ellip_acd_cmplx ; +p_ellip_ap_zp ellip_ap_zp ; p_ellip_asn ellip_asn ; p_ellip_asn_cmplx ellip_asn_cmplx ; p_ellip_cd ellip_cd ; @@ -177,6 +178,7 @@ void* dspl_load() LOAD_FUNC(dmod); LOAD_FUNC(ellip_acd); LOAD_FUNC(ellip_acd_cmplx); + LOAD_FUNC(ellip_ap_zp); LOAD_FUNC(ellip_asn); LOAD_FUNC(ellip_asn_cmplx); LOAD_FUNC(ellip_cd); diff --git a/include/dspl.h b/include/dspl.h index 19d7e62..979652d 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -310,6 +310,14 @@ DECLARE_FUNC(int, ellip_acd_cmplx, complex_t* w COMMA double k COMMA complex_t* u); //------------------------------------------------------------------------------ +DECLARE_FUNC(int, ellip_ap_zp, int ord + COMMA double rp + COMMA double rs + COMMA complex_t* z + COMMA int* nz + COMMA complex_t* p + COMMA int* np); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, ellip_asn, double* w COMMA int n COMMA double k diff --git a/release/include/dspl.c b/release/include/dspl.c index a3fc94e..c4e57e4 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -56,6 +56,7 @@ p_dft_cmplx dft_cmplx ; p_dmod dmod ; p_ellip_acd ellip_acd ; p_ellip_acd_cmplx ellip_acd_cmplx ; +p_ellip_ap_zp ellip_ap_zp ; p_ellip_asn ellip_asn ; p_ellip_asn_cmplx ellip_asn_cmplx ; p_ellip_cd ellip_cd ; @@ -177,6 +178,7 @@ void* dspl_load() LOAD_FUNC(dmod); LOAD_FUNC(ellip_acd); LOAD_FUNC(ellip_acd_cmplx); + LOAD_FUNC(ellip_ap_zp); LOAD_FUNC(ellip_asn); LOAD_FUNC(ellip_asn_cmplx); LOAD_FUNC(ellip_cd); diff --git a/release/include/dspl.h b/release/include/dspl.h index 19d7e62..979652d 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -310,6 +310,14 @@ DECLARE_FUNC(int, ellip_acd_cmplx, complex_t* w COMMA double k COMMA complex_t* u); //------------------------------------------------------------------------------ +DECLARE_FUNC(int, ellip_ap_zp, int ord + COMMA double rp + COMMA double rs + COMMA complex_t* z + COMMA int* nz + COMMA complex_t* p + COMMA int* np); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, ellip_asn, double* w COMMA int n COMMA double k