added functions for butterworth filter

dspl/src/filter.c

@@ -1,71 +0,0 @@
-/*
-* 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;
-    
-}
-
-
-
-
-
-

dspl/src/filter_ap.c

@@ -0,0 +1,168 @@
+/*
+* 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 <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "dspl.h"
+
+
+
+
+
+/**************************************************************************************************
+Analog Normalized Butterworth filter  
+***************************************************************************************************/
+int DSPL_API butter_ap(double Rp, int ord, double* b, double* a)
+{
+	int res;
+    complex_t *z = NULL;
+    complex_t *p = NULL;
+    int nz, np;
+
+    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 = butter_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;
+
+    b[0] = a[0];
+
+
+exit_label:
+    if(z)
+        free(z);
+    if(p)
+        free(p);
+    return res;   
+}
+
+
+
+
+
+/**************************************************************************************************
+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 ep; 
+    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);
+
+    alpha = pow(ep, -1.0/(double)ord);
+    if(r)
+    {
+        RE(p[ind]) = -alpha;
+        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]) = -alpha *  sin(theta);
+        IM(p[ind]) =   alpha *  cos(theta);
+        IM(p[ind+1]) =  -alpha *  cos(theta);
+        ind+=2;
+    }  
+    *np = ord;
+    *nz = 0;  
+    return RES_OK;
+}
+
+
+
+
+
+
+
+
+/**************************************************************************************************
+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;
+    
+}
+
+
+
+
+
+

dspl/src/polyval.c

@@ -35,7 +35,7 @@ int DSPL_API poly_z2a_cmplx(complex_t* z, int nz, int ord, complex_t* a)
 
     if(!z || !a)
         return ERROR_PTR;
-    if(nz < 1)
+    if(nz < 0)
         return ERROR_SIZE;
     if(nz > ord || ord < 1)
         return ERROR_POLY_ORD;

include/dspl.c

@@ -35,10 +35,10 @@
 
 #ifndef BUILD_LIB
 
-
-
 p_acos_cmplx                acos_cmplx          ;
 p_asin_cmplx                asin_cmplx          ;
+p_butter_ap                 butter_ap           ;                
+p_butter_ap_zp              butter_ap_zp        ;
 p_cheby_poly1               cheby_poly1         ;
 p_cheby_poly2               cheby_poly2         ;
 p_cmplx2re                  cmplx2re            ;
@@ -135,6 +135,8 @@ void* dspl_load()
 
     LOAD_FUNC(acos_cmplx);
     LOAD_FUNC(asin_cmplx);
+    LOAD_FUNC(butter_ap);
+    LOAD_FUNC(butter_ap_zp);
     LOAD_FUNC(cheby_poly1);
     LOAD_FUNC(cheby_poly2);
     LOAD_FUNC(cmplx2re);

include/dspl.h

@@ -194,6 +194,10 @@ extern "C" {
 
 DECLARE_FUNC(int,   acos_cmplx,         complex_t*  COMMA int       COMMA complex_t*);
 DECLARE_FUNC(int,   asin_cmplx,         complex_t*  COMMA int       COMMA complex_t*);
+
+DECLARE_FUNC(int,   butter_ap,          double      COMMA int       COMMA double* COMMA double*);
+DECLARE_FUNC(int,   butter_ap_zp,       int         COMMA double    COMMA complex_t* COMMA int* COMMA complex_t* COMMA int*);
+
 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,   cmplx2re,           complex_t*  COMMA int       COMMA double* COMMA double*);

release/include/dspl.c

@@ -35,10 +35,10 @@
 
 #ifndef BUILD_LIB
 
-
-
 p_acos_cmplx                acos_cmplx          ;
 p_asin_cmplx                asin_cmplx          ;
+p_butter_ap                 butter_ap           ;                
+p_butter_ap_zp              butter_ap_zp        ;
 p_cheby_poly1               cheby_poly1         ;
 p_cheby_poly2               cheby_poly2         ;
 p_cmplx2re                  cmplx2re            ;
@@ -135,6 +135,8 @@ void* dspl_load()
 
     LOAD_FUNC(acos_cmplx);
     LOAD_FUNC(asin_cmplx);
+    LOAD_FUNC(butter_ap);
+    LOAD_FUNC(butter_ap_zp);
     LOAD_FUNC(cheby_poly1);
     LOAD_FUNC(cheby_poly2);
     LOAD_FUNC(cmplx2re);

release/include/dspl.h

@@ -194,6 +194,10 @@ extern "C" {
 
 DECLARE_FUNC(int,   acos_cmplx,         complex_t*  COMMA int       COMMA complex_t*);
 DECLARE_FUNC(int,   asin_cmplx,         complex_t*  COMMA int       COMMA complex_t*);
+
+DECLARE_FUNC(int,   butter_ap,          double      COMMA int       COMMA double* COMMA double*);
+DECLARE_FUNC(int,   butter_ap_zp,       int         COMMA double    COMMA complex_t* COMMA int* COMMA complex_t* COMMA int*);
+
 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,   cmplx2re,           complex_t*  COMMA int       COMMA double* COMMA double*);

test/bin/gnuplot/butter_ap_test.plt

@@ -0,0 +1,7 @@
+
+set logscale x
+unset key
+set grid
+set xlabel "frequency, rad/s"
+set ylabel "Butterworth filter magnitude, dB"
+plot 'dat/butter_ap_test_mag.txt' with lines

test/src/butter_ap_test.c

@@ -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 = butter_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/butter_ap_test_mag.txt");    
+    
+    dspl_free(handle);      // free dspl handle
+
+    system("gnuplot -p  gnuplot/butter_ap_test.plt");
+
+    return 0;
+}
+
+