gr-filter: Remove simpson's rule integration and replace with GSL call

author: Andy Walls <awalls.cx18@gmail.com> 2018-02-08 19:58:00 -0500
committer: Marcus Müller <marcus@hostalia.de> 2018-02-17 20:08:44 +0100
commit: caa4e5d34daf8e27c6ae7ddf982ae25113a934e3 (patch)
tree: 41920df50ec74bc6e6749689edb58dc124dd3c0f
parent: 1e5bd96c17c5c454d7a75450c2d54f18d12392b0 (diff)
5 files changed, 29 insertions, 86 deletions
diff --git a/gr-filter/lib/gen_interpolator_taps/gen_interpolator_taps.c b/gr-filter/lib/gen_interpolator_taps/gen_interpolator_taps.c
index 2f359102cf..484899e6e5 100644
--- a/gr-filter/lib/gen_interpolator_taps/gen_interpolator_taps.c
+++ b/gr-filter/lib/gen_interpolator_taps/gen_interpolator_taps.c
@@ -1,6 +1,6 @@
 /* -*- c++ -*- */
 /*
- * Copyright 2002 Free Software Foundation, Inc.
+ * Copyright 2002-2018 Free Software Foundation, Inc.
  *
  * This file is part of GNU Radio
  *
@@ -23,6 +23,7 @@
 #include <stdio.h>
 #include <unistd.h>
 #include <stdlib.h>
+#include <gsl/gsl_integration.h>
 
 #define	NSTEPS		 10	// how many steps of mu are in the generated table
 #define	MAX_NSTEPS	256
@@ -33,6 +34,7 @@ extern void initpt (double x[], int ntaps);
 extern double objective (double x[], int *ntaps);
 extern double global_mu;
 extern double global_B;
+extern gsl_integration_workspace *global_gsl_int_workspace;
 
 // fortran
 extern double prax2_ (double (fct)(double x[], int *ntaps),
@@ -123,6 +125,13 @@ main (int argc, char **argv)
     exit (1);
   }
 
+  global_gsl_int_workspace = gsl_integration_workspace_alloc(4000);
+  if (global_gsl_int_workspace == NULL) {
+    fprintf (stderr, "%s: unable to allocate GSL integration work space\n",
+             argv[0]);
+    exit (1);
+  }
+
   step_size = 1.0/nsteps;
 
   // the optimizer chokes on the two easy cases (0/N and N/N).  We do them by hand...
@@ -152,6 +161,8 @@ main (int argc, char **argv)
     }
   }
 
+  gsl_integration_workspace_free(global_gsl_int_workspace);
+
   // now compute remaining values via symmetry
 
   for (j = 0; j < nsteps/2; j++){
diff --git a/gr-filter/lib/gen_interpolator_taps/generate.sh b/gr-filter/lib/gen_interpolator_taps/generate.sh
index a7a6c6cafe..82ade73cee 100755
--- a/gr-filter/lib/gen_interpolator_taps/generate.sh
+++ b/gr-filter/lib/gen_interpolator_taps/generate.sh
@@ -31,12 +31,12 @@ gfortran -o gen_interp_differentiator_taps \
 ./gen_interp_differentiator_taps -n 128 -t 8 -B 0.25 \
         > interp_differentiator_taps.h
 
-gcc      -c -o simpson.o               simpson.c
 gcc      -c -o objective_fct.o         objective_fct.c
 gcc      -c -o gen_interpolator_taps.o gen_interpolator_taps.c
 
 gfortran -o gen_interpolator_taps \
-        gen_interpolator_taps.o praxis.o objective_fct.o simpson.o
+        gen_interpolator_taps.o praxis.o objective_fct.o \
+        -lgsl -lgslcblas
 
 ./gen_interpolator_taps -n 128 -t 8 -B 0.25 \
         > interpolator_taps.h
diff --git a/gr-filter/lib/gen_interpolator_taps/objective_fct.c b/gr-filter/lib/gen_interpolator_taps/objective_fct.c
index 129486d634..061fdf4970 100644
--- a/gr-filter/lib/gen_interpolator_taps/objective_fct.c
+++ b/gr-filter/lib/gen_interpolator_taps/objective_fct.c
@@ -1,6 +1,6 @@
 /* -*- c -*- */
 /*
- * Copyright 2002 Free Software Foundation, Inc.
+ * Copyright 2002-2018 Free Software Foundation, Inc.
  *
  * This file is part of GNU Radio
  *
@@ -26,7 +26,7 @@
 
 #include <math.h>
 #include <assert.h>
-#include "simpson.h"
+#include <gsl/gsl_integration.h>
 
 #define	MU	0.5			/* the MU for which we're computing coeffs */
 
@@ -39,6 +39,8 @@ static double  *global_h;
 double		global_mu = MU;
 double		global_B  = B;
 
+gsl_integration_workspace *global_gsl_int_workspace = NULL;
+
 /*
  * This function computes the difference squared between the ideal
  * interpolator frequency response at frequency OMEGA and the
@@ -49,7 +51,7 @@ double		global_B  = B;
  */
 
 static double
-integrand (double omega)
+integrand (double omega, void *params)
 {
   double real_ideal;
   double real_approx;
@@ -89,10 +91,19 @@ integrand (double omega)
 double
 c_fcn (double *x, int n)
 {
+  gsl_function F;
+  double result, error;
+
+  F.function = integrand;
+  F.params = NULL;
+
   assert ((n & 1) == 0);	/* assert n is even */
   global_n = n;
   global_h = x;
-  return qsimp (integrand, -2 * M_PI * global_B, 2 * M_PI * global_B);
+  gsl_integration_qag(&F, -2 * M_PI * global_B, 2 * M_PI * global_B,
+                      0.0, 1e-12, 1000, GSL_INTEG_GAUSS61,
+                      global_gsl_int_workspace, &result, &error);
+  return result;
 }
 
 /* this is the interface expected by the calling fortran code */
diff --git a/gr-filter/lib/gen_interpolator_taps/simpson.c b/gr-filter/lib/gen_interpolator_taps/simpson.c
deleted file mode 100644
index 31aaae4aef..0000000000
--- a/gr-filter/lib/gen_interpolator_taps/simpson.c
+++ /dev/null
@@ -1,76 +0,0 @@
-/* -*- c -*- */
-#include <math.h>
-#include <stdio.h>
-
-#define	EPS	(1.0e-5)
-#define	JMAX	16
-
-/*
- * Compute the Nth stage of refinement of an extended trapezoidal
- * rule.  FUNC is input as a pointer to a function to be integrated
- * between limits A and B.  When called with N = 1, the routine
- * returns the crudest estimate of the integral from A to B of f(x)
- * dx.  Subsequent calls with N=2,3,... (in that sequential order)
- * will improve the accuracy by adding 2**(N-2) additional interior
- * points.
- *
- * N.B., this function contains static state and IS NEITHER RENTRANT
- * NOR THREAD SAFE!
- */
-
-double
-trapzd (double (*func)(double),
-	double a, double b,
-	int    n)
-{
-  long double		x, tnm, sum, del;
-  static long double	s;
-  static int	it;
-  int		j;
-
-  if (n == 1){
-    it = 1;	/* # of points to add on the next call */
-    s = 0.5 * (b - a) * (func(a) + func(b));
-    return s;
-  }
-  else {
-    tnm = it;
-    del = (b-a)/tnm;		/* this is the spacing of the points to be added */
-    x = a + 0.5*del;
-    for (sum = 0.0, j = 1; j <= it; j++, x += del)
-      sum += func(x);
-    it *= 2;
-    s = 0.5 * (s + (b-a) * sum/tnm);	/* replace s by it's refined value */
-    return s;
-  }
-}
-
-/*
- * Returns the integral of the function FUNC from A to B.  The
- * parameters EPS can be set to the desired fractional accuracy and
- * JMAX so that 2**(JMAX-1) is the maximum allowed number of steps.
- * Integration is performed by Simpson's rule.
- */
-
-double
-qsimp (double (*func)(double),
-       double a,	/* lower limit */
-       double b)	/* upper limit */
-{
-  int		j;
-  long double	s, st, ost, os;
-
-  ost = os = -1.0e30;
-  for (j = 1; j <= JMAX; j++){
-    st = trapzd (func, a, b, j);
-    s = (4.0 * st - ost)/3.0;
-    if (fabs (s - os) < EPS * fabs(os))
-      return s;
-    os = s;
-    ost = st;
-  }
-  fprintf (stderr, "Too many steps in routine QSIMP\n");
-  // exit (1);
-  return s;
-}
-
diff --git a/gr-filter/lib/gen_interpolator_taps/simpson.h b/gr-filter/lib/gen_interpolator_taps/simpson.h
deleted file mode 100644
index 68774f9a2e..0000000000
--- a/gr-filter/lib/gen_interpolator_taps/simpson.h
+++ /dev/null
@@ -1,3 +0,0 @@
-double qsimp (double (*func)(double),
-	      double a, double b);
-
author	Andy Walls <awalls.cx18@gmail.com>	2018-02-08 19:58:00 -0500
committer	Marcus Müller <marcus@hostalia.de>	2018-02-17 20:08:44 +0100
commit	caa4e5d34daf8e27c6ae7ddf982ae25113a934e3 (patch)
tree	41920df50ec74bc6e6749689edb58dc124dd3c0f
parent	1e5bd96c17c5c454d7a75450c2d54f18d12392b0 (diff)