From 6fb251d989654388dd32166e7cfd982306a882fa Mon Sep 17 00:00:00 2001
From: Dsplib <admin@dsplib.org>
Date: Wed, 28 Mar 2018 00:02:37 +0300
Subject: [PATCH] added functions poly_z2a and poly_zp2ab

---
 dspl/src/filter.c      |  71 +++++++++++++++++++++++++++++++++++
 dspl/src/polyval.c     |  83 ++++++++++++++++++-----------------------
 include/dspl.c         |   4 ++
 include/dspl.h         |   3 +-
 release/include/dspl.c |   4 ++
 release/include/dspl.h |   3 +-
 test/bin/dft_test.exe  | Bin 14512 -> 14584 bytes
 7 files changed, 119 insertions(+), 49 deletions(-)
 create mode 100644 dspl/src/filter.c

diff --git a/dspl/src/filter.c b/dspl/src/filter.c
new file mode 100644
index 0000000..7dc514b
--- /dev/null
+++ b/dspl/src/filter.c
@@ -0,0 +1,71 @@
+/*
+* 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 "dspl.h"
+
+
+
+
+
+/**************************************************************************************************
+Zeros and poles to filter coefficients recalc
+***************************************************************************************************/
+int DSPL_API filter_zp2ab(complex_t *z, int nz, complex_t *p, int np, int ord, double* b, double* a)
+{
+    complex_t *acc = NULL;
+    int res;
+
+    if(!z || !p || !b || !a)
+        return ERROR_PTR;
+    if(nz < 0 || np < 0)
+        return ERROR_SIZE;
+    if(nz < ord || np < ord)
+        return ERROR_POLY_ORD;
+
+    acc = (complex_t*) malloc((ord+1) * sizeof(complex_t));
+    res = poly_z2a_cmplx(z, nz, ord, acc);
+    if(res != RES_OK)
+        goto exit_label;    
+    res = cmplx2re(acc, ord+1, b, NULL);
+    if(res != RES_OK)
+        goto exit_label;
+
+    res = poly_z2a_cmplx(p, np, ord, acc);
+    if(res != RES_OK)
+        goto exit_label;    
+    res = cmplx2re(acc, ord+1, a, NULL);
+    if(res != RES_OK)
+        goto exit_label;
+
+exit_label:
+    if(acc)
+        free(acc);
+    return res;
+    
+}
+
+
+
+
+
+
diff --git a/dspl/src/polyval.c b/dspl/src/polyval.c
index d80cc66..10a184b 100644
--- a/dspl/src/polyval.c
+++ b/dspl/src/polyval.c
@@ -25,6 +25,42 @@
 
 
 
+/**************************************************************************************************
+Polynomial zeros to coefficients
+***************************************************************************************************/
+int DSPL_API poly_z2a_cmplx(complex_t* z, int nz, int ord, complex_t* a)
+{
+    int k, ind, res;
+    complex_t x[2];
+
+    if(!z || !a)
+        return ERROR_PTR;
+    if(nz < 1)
+        return ERROR_SIZE;
+    if(nz > ord || ord < 1)
+        return ERROR_POLY_ORD;
+
+    RE(x[1]) = 1.0;
+    IM(x[1]) = 0.0;
+
+    memset(a, 0, (ord+1) * sizeof(complex_t));
+
+    RE(a[0]) = 1.0;
+    ind = 1;
+    for(k = 0; k < nz; k++)
+    {
+        RE(x[0]) = -RE(z[k]);
+        IM(x[0]) = -IM(z[k]);
+        res = conv_cmplx(a, ind, x, 2, a);
+        if(res!=RES_OK)
+            return res;
+        ind++;
+    }
+
+    return RES_OK;
+}
+
+
 
 
 /**************************************************************************************************
@@ -88,50 +124,3 @@ int DSPL_API polyval_cmplx(complex_t* a, int ord, complex_t* x, int n, complex_t
 
 
 
-
-
-/**************************************************************************************************
-polynomial zeros to coeff reecalc 
-***************************************************************************************************/
-int poly_z2a(complex_t *z, int nz, int ord, complex_t *a)
-{
-    int k, ind, res;
-    complex_t x[2];
-
-    if(!z || !a)
-	    return ERROR_PTR;	
-	if(ord<0)
-		return ERROR_POLY_ORD;	
-	if(nz<1 || nz > ord)
-		return ERROR_SIZE;
-
-    memset(a, 0, (ord+1) * sizeof(complex_t));
-    RE(a[0]) = 1.0;
-    IM(a[0]) = 0.0;
-
-    RE(x[1]) = 1.0;
-    IM(x[1]) = 0.0;
-
-
-    ind = 1;
-    for(k=0; k < nz; k++)
-    {
-        RE(x[0]) = -RE(z[k]);
-        IM(x[0]) = -IM(z[k]);
-        res = conv_cmplx(a, ind, x, 2, a);
-        if(res!=RES_OK)
-            return res;        
-        ind++;        
-    }
-    return RES_OK;
-}
-
-
-
-
-
-
-
-
-
-
diff --git a/include/dspl.c b/include/dspl.c
index ec12f5b..d81c364 100644
--- a/include/dspl.c
+++ b/include/dspl.c
@@ -52,6 +52,7 @@ p_dmod                      dmod                ;
 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          ; 
@@ -66,6 +67,7 @@ p_goertzel_cmplx            goertzel_cmplx      ;
 p_linspace                  linspace            ;
 p_log_cmplx                 log_cmplx           ;
 p_logspace                  logspace            ;
+p_poly_z2a_cmplx            poly_z2a_cmplx      ;
 p_polyval                   polyval             ;
 p_polyval_cmplx             polyval_cmplx       ;
 p_randn                     randn               ;
@@ -144,6 +146,7 @@ void* dspl_load()
     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);
@@ -158,6 +161,7 @@ void* dspl_load()
     LOAD_FUNC(linspace);
     LOAD_FUNC(log_cmplx);
     LOAD_FUNC(logspace);
+    LOAD_FUNC(poly_z2a_cmplx);
     LOAD_FUNC(polyval);
     LOAD_FUNC(polyval_cmplx);
     LOAD_FUNC(randn);
diff --git a/include/dspl.h b/include/dspl.h
index 8805244..4513d74 100644
--- a/include/dspl.h
+++ b/include/dspl.h
@@ -207,7 +207,7 @@ DECLARE_FUNC(double,dmod,               double      COMMA double);
 DECLARE_FUNC(int,   farrow_lagrange,    double*     COMMA int       COMMA double COMMA double COMMA double COMMA double** COMMA int*);
 DECLARE_FUNC(int,   farrow_spline,      double*     COMMA int       COMMA double COMMA double COMMA 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* COMMA complex_t* );
 DECLARE_FUNC(int,   fft_cmplx,          complex_t*  COMMA int       COMMA fft_t* COMMA complex_t* );
@@ -223,6 +223,7 @@ DECLARE_FUNC(int,   goertzel_cmplx,     complex_t*  COMMA int       COMMA int* C
 DECLARE_FUNC(int,   linspace,           double      COMMA double    COMMA int COMMA int COMMA double*);
 DECLARE_FUNC(int,   log_cmplx,          complex_t*  COMMA int       COMMA complex_t*);
 DECLARE_FUNC(int,   logspace,           double      COMMA double    COMMA int COMMA int COMMA double*);
+DECLARE_FUNC(int,   poly_z2a_cmplx,     complex_t*  COMMA int       COMMA int COMMA complex_t*);
 DECLARE_FUNC(int,   polyval,            double*     COMMA int       COMMA double* COMMA int COMMA double*);
 DECLARE_FUNC(int,   polyval_cmplx,      complex_t*  COMMA int       COMMA complex_t* COMMA int COMMA complex_t*);
 DECLARE_FUNC(int,   randn,              double*     COMMA int       COMMA double COMMA double);
diff --git a/release/include/dspl.c b/release/include/dspl.c
index ec12f5b..d81c364 100644
--- a/release/include/dspl.c
+++ b/release/include/dspl.c
@@ -52,6 +52,7 @@ p_dmod                      dmod                ;
 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          ; 
@@ -66,6 +67,7 @@ p_goertzel_cmplx            goertzel_cmplx      ;
 p_linspace                  linspace            ;
 p_log_cmplx                 log_cmplx           ;
 p_logspace                  logspace            ;
+p_poly_z2a_cmplx            poly_z2a_cmplx      ;
 p_polyval                   polyval             ;
 p_polyval_cmplx             polyval_cmplx       ;
 p_randn                     randn               ;
@@ -144,6 +146,7 @@ void* dspl_load()
     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);
@@ -158,6 +161,7 @@ void* dspl_load()
     LOAD_FUNC(linspace);
     LOAD_FUNC(log_cmplx);
     LOAD_FUNC(logspace);
+    LOAD_FUNC(poly_z2a_cmplx);
     LOAD_FUNC(polyval);
     LOAD_FUNC(polyval_cmplx);
     LOAD_FUNC(randn);
diff --git a/release/include/dspl.h b/release/include/dspl.h
index 8805244..4513d74 100644
--- a/release/include/dspl.h
+++ b/release/include/dspl.h
@@ -207,7 +207,7 @@ DECLARE_FUNC(double,dmod,               double      COMMA double);
 DECLARE_FUNC(int,   farrow_lagrange,    double*     COMMA int       COMMA double COMMA double COMMA double COMMA double** COMMA int*);
 DECLARE_FUNC(int,   farrow_spline,      double*     COMMA int       COMMA double COMMA double COMMA 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* COMMA complex_t* );
 DECLARE_FUNC(int,   fft_cmplx,          complex_t*  COMMA int       COMMA fft_t* COMMA complex_t* );
@@ -223,6 +223,7 @@ DECLARE_FUNC(int,   goertzel_cmplx,     complex_t*  COMMA int       COMMA int* C
 DECLARE_FUNC(int,   linspace,           double      COMMA double    COMMA int COMMA int COMMA double*);
 DECLARE_FUNC(int,   log_cmplx,          complex_t*  COMMA int       COMMA complex_t*);
 DECLARE_FUNC(int,   logspace,           double      COMMA double    COMMA int COMMA int COMMA double*);
+DECLARE_FUNC(int,   poly_z2a_cmplx,     complex_t*  COMMA int       COMMA int COMMA complex_t*);
 DECLARE_FUNC(int,   polyval,            double*     COMMA int       COMMA double* COMMA int COMMA double*);
 DECLARE_FUNC(int,   polyval_cmplx,      complex_t*  COMMA int       COMMA complex_t* COMMA int COMMA complex_t*);
 DECLARE_FUNC(int,   randn,              double*     COMMA int       COMMA double COMMA double);
diff --git a/test/bin/dft_test.exe b/test/bin/dft_test.exe
index e4c891e8a7549283084833468a5a7a0fce9d711d..7f29d76dc4cd2f1f02fd2921a9dfa1720d489c85 100755
GIT binary patch
delta 2748
zcmZA33rv$&6bJD8wLID)+<q-kL6o9JDJW2bf-_U<P>M4SMJK^!R7K(gwGf;RU5e8i
z?KUFy>gLRJa~gGSnZwSEp{W?%1QWL?(M(Myi`m3YH{5*Uqr0c~=rTy?ef)pt^!EE&
z%2%}yYMV1O1rC4aKK<5s68_ohzm%_Za8=o8B4)Wo%!P7d<_DF_@gI~Qel&C2;Fj-`
z#_e13_PEYbht7Kf+cmxj?@~P>CB9RxrmeO{Q4utX2bb?UwMW>hX-4F0BBu5tn-Lb3
zxUj-vJ}z=$u@D!9uvmnPQdmHLXyAs$a$I;q2}_f#gr<@R-|xYCzBhWLW`K^Lk1o`?
zjfKQ@`a~~J*4;{IGZm7a#wVeWi?!1Y`=j7|Do@c5)*018QlA$5g*)|D?;x>o`L6bq
zf%HYLwOC2d99h){*~(+V;Ye~NSqCH!31llj430(;E6JZgR3wm5Ru>9|f)kMa_Bi}x
z#TF~?)*JF*o3fWd6r(sUMR%fTbQ0#Qy$v$TY`~_ge25`AZ+dv=$y2dYhK04V@E{bf
z#6rj@y9i8Gd5t02tA=+D!OrVgI7}8+LE#82gp9IHz~&=+{3v|%v%m;4$S9i!tOnTy
z$#wy|1aonI$S8Y3FmJFH-IvnMNf7QZI>;z9L%VgzHcIA!cJCsCjIx{JDhD?sb4xY@
zx>X^AjI!^5?LaopE8RhGFQS8xvR&ZzqccgSg=bAf1{r0mfPIeakG}9F`3SoG1#|IQ
zA)~Am*b!u>B)brACElJ6bP!Tz2iJ~nr*w(XZZ|T>DANGzMD~VcBSgEN24s*?c17Hl
zU>CAiB+~&aL)PUDyRX4rMwcO7C%BR5Afv1a*mY!zWK)4fB7=;wYG8fHE<X%kk_Eti
zhPk*TkWn@fm^W5@IN?(lknRg`-=KqxvVq_dk-aBbEU;!|kWqG9+?HTEvW1cb&4MjP
z1{q~%fMp|_@X#AR-3oA%&_PJq9&iQdEV3N|HV7GHl&uBkMArKt-0d@9cVI4FD`b?-
z23Cr!Q?d$R-m~Z+q%0d;Il51!+hY=LA2P@&iw0JWY>i~La64p>QPwSPORxsnEXi&`
zw{m1PCNH>e!PTJ~CEXNo4s?)FwguQ4WKoh)U<PE6QMMFViYjmG)%)Q~@(HkOFc+5u
zGRlg8ZS-QhR_XlUj-lHaZZ{a*7G!?OW~zd1Lk1aT_rz@p?m)Igvhl!{BZG{xAAub}
z_R{^37yi@_09T9-Ldp(+YeAPH+a&<AN>-oL2I-AP;hK-17So>h?}b8+_PHlrzMtlv
z6kkf!*)(;7_#PzW#!(Udw`_LKalY#Gg`BgUvxZG;N`D$_C9_>kqhKs6@!c!&U6^M7
z)wTYn&UNC6wkJC91AR-0SzJ@Uu&%PfGkV^9UZQ?%sB9QLuX4dMPt9t6O+7NHNIXxp
z_$~cM97v44Vr&&7xf=?(VH8K`6+^f9G!0ZlF2-Cjy2W3T)nX)sUo;KXY`??%Oxc>&
zzj>y4xRT3fnDgyA(V|OyYZ7OialJ-UVTjXBj?pfQ5~DMwC{7nIqP^m`owt}RN-OU)
z7e{uBcc><i7QaTnWk{(~6__#PxK^p+>QIZ~5swKeDtZ3Tg(!XYqa`NdG6yP##%Yvh
zzSy=|L(II-maVu0_Vj2)G;$258qva;k)!nTTJhMycW10nqIs>|q8Pbqw<M)&2{Dpn
zvdHRw);BylxKo(1VN)uMDvT>l8Xd|?zFTyt!?e6#+Q84*bCf3T$h0JJVN4{Bbko9c
zZ6tQs$^DrYM+<gX-LJz@+;8l+e}<Q2rmCCz?Ox<3GA+s#=$LdDE4TMEA%0sJ15eAc
zD2aSsmL<{_Ax_56DQ0`Vn8B-T?N>7%gg86V%d>a!lC0Fo$<ONa6<F6E@MRga|9k8C
zvybRJD{9_^w_jjl&MhU7>u_t8I_}TQ(N`H|zbf9Dm#=zo<uS@j)5O&=Ifm!A-OQ(s
t8EbIiznJ2#zcVIEqqw<0f0ZUHmm3PEXfk5>%!0A{jx5<=0rwZA{tKy>RcQbK

delta 2612
zcmZA3drVVT90%}wZ+W!U^|Yl}1*L^mrGOwAh%TG98*14M!`E;=GK`PJ$L83;$A+n!
zo@rG!;Nap+TyQhS_=h8_*^sdU!Da3}O*Ej1k~tl2%xQchdiVR?^K*Z2liqvJ`F?)C
z`|9gnWpia+rgnW|ebyn1&!rfPufC^9Z4nN2yR*?Lry89BW2ez6j;QBcgO|SCbEb7x
z)$L~6l9r}hEwfF_6QX;yNL)D8q9`Si6QSBh-xjEd*<lrm{K@=rzLhL9emKO;U+lBv
zfYHE<11Ak8;vks@0UV5^K^YFPKQ^eqK^hI-(397U@+!=u#6|v!`NYA*9IY);oK6gw
z?e>5Y`mVD{q?rFnE`_V;+wxG?L%f4+do=;hufr}zTgI3>fK{c(szuP^4gZKWk&AtE
zQAXlx++L-xTvoMWR^1jmOJt%X?;<H6f?4&)vEPY2A;|?K4kDN_{-&<$v1`Px3>pWU
z>Gq0Vi?sl!8GjMQT@vk5^p1c=UAWIXS}|kXhwPDEq*+r69vL$aAHj^tA+MuS79PRE
zYAVEx@w3P}?P9GpCERHgevUIXDtuBFF2}<0REQblJCHev4L&eV{TpNoG0YfWiflBo
zGm?FPtOxhu{4rzvX<*@4I=Q3Ll_!Jyf*fXyCt<sMV%sIFz;?TcVa9kr+~wFbVil6j
z#%{}rVaE8k$X*~eB`jSEwtJQwhKzrLu9Tcpvahk-2x6Eqz5&?+Vt)@9w`9-->|fl6
zYsHN5GGxn$osg^$x`q2rlEaX3Kf1N#_Dffa?G6ybjByRwdSdG(%R}}yG0Yf05053b
znb=E`^&y)}Y;)Lf$I$I0mnmH(x*T$tF<y&oKQUFZ_mIUC!;JBj$eM|ryKmf*E@T&Q
zA8rX|j88)rK29zwT@Jdh$YI9#NOUKO?Uw9v60kaAm@$3>9!u;Lu@#biglrWt%ozUx
z*?D49?}v@E8;|Zeau_oH8M;5nxn;X1WU0h3W4s1gKe6lgjBYoaz;59_Tq|ab&qp>u
z?0d;pW4G|n<S=AB7o7s%SNKKrsdP`FJ46gK#uJex5PL_mMm%bm7-o$3!efasV)G<h
zf@}dX<_x1dj_xPBJn97LE}+XNhZ*Ay$Wp0Yf@HrUvl7FM@ioZC5xaQTxFt;vV3%+o
zZV6_L7bD9JlWUPKg078RuHnX@dxluOWDAk)A%+>_ci^$arW0E&*)n9WOIGD+#q{bg
z;Ft%0q?GnW&cGaaSKiU}{EkIkp-A_lF6Jp#l&iy_W@Rq(Y(u8h1<Qjk2P1lLesJEb
z+OefJuQESWI{~+IOConlB4=j#FNe1Fn?s!sO?~CjH2Y>%9Cqv%^PJ671JEE5UTrsZ
zeK%|u-_rGCuw8pw*B9d!j@}R3L3o$DR7DvMpX?mi*1)z(6*H1DH1ln-Iw@Dnyeal2
zjZ>RMXHt=WBb4ugUos_Vvv1Y3g;tw+`Y_Yj1lR_L72C`%hz`Q5SERe#YPXo;ni-E@
z(pt3Mway}%Tr~}8WvVMWJFVTMx`ZRct$JadJldf={?Ui({}+?uRJxtfg&8&t+$!IP
z8vG>~@a3xZsDEst3XSr|F^qzld2oi6ux=3tGS{glQR#Q9BSp~f9_2GB(8i`*PY4)|
zt#FyR_zm;CaDbTu#-XGrq4WwPGaG9~Ss@NU`WmGRhNN$aQ~tbB?^FLe7BH&fskTn!
zXSvmbLWT8Fu_?>r*f^xVRV;<99bHKe9b?B3bBX5bhRk21`MpEt2gHpmcf2_cu34N?
zll-Z$hgWuZi1VO;yP{a+S?>~?vOV$s$93vdk)G|2p3HuriQjUg(T2P$s_0eYq8sxo
zOlpScDaf<9;Nyu4JjcQ!hnGfPyD_O$idPErtdE~6Jhz6zN!D)oEZ_nU^>kr^rk07G
fq7B-`JTYhT4DFs(?3+Bvx;NV(%FH;iFF)-+jrbu?