/*----------------------------------------------------------------------
bfilt - High-order Butterworth filter.
by Andy Allinger, 2021, released to the public domain

   Permission  to  use, copy, modify, and distribute this software and
   its documentation  for  any  purpose  and  without  fee  is  hereby
   granted,  without any conditions or restrictions.  This software is
   provided "as is" without express or implied warranty.

Refer to:

    "Cookbook formulae for audio EQ biquad filter coefficients"
    Robert Bristow-Johnson, [2005]

    "The Butterworth Low-Pass Filter"
    John Stensby, 19 Oct 2005

The Butterworth poles lie on a circle.  The product of the qualities
is 2^(-1/2).

___Name____Type______In/Out___Description_______________________________
   s[l]    double*   Both     Sample data, returned filtered
   l       int       In       Length of s
   fsamp   double    In       Sample rate, cycles per second
   fpass   double    In       Pass frequency, cycles per second
   fstop   double    In       Stop frequency, cycles per second
   hpass   double    In       Minimum passband transmission, fraction 0...1
   hstop   double    In       Maximum stopband transmission, fraction 0...1
----------------------------------------------------------------------*/
#include <stdbool.h>
#include <math.h>

int bfilt (double *s, int l, double fsamp,
           double fpass, double fstop, double hpass, double hstop)
{
/* Local variables */
	int j, k, n;
	double a0, a1, a2, b1, b2, c, d, e, fcut, o, q, r, w0, x, x1, x2, y, y1, y2;
	bool lo;
/* Function Body */
	if (l <= 0) return 1;
	if (fsamp <= 0.) return 2;
	if (fpass <= 0. || fpass >= 0.5 * fsamp) return 3;
	if (fstop <= 0. || fstop >= 0.5 * fsamp) return 4;
	if (hpass <= 0. || hpass >= 1.) return 5;
	if (hstop <= 0. || hstop >= 1.) return 6;
	if (fpass == fstop) return 7;
/* Find coefficients */
	lo = fpass < fstop;
	d = 1. / hstop;
	e = sqrt(1. / (hpass * hpass) - 1.);
	n = ((int) fabs(log(e / sqrt(d * d - 1.)) / log(fpass / fstop))) + 1;
	if (n % 2) ++n;
	o = lo ? -1. / (double) n : +1. / (double) n;
	fcut = fstop * pow(sqrt(d * d - 1.), o);
	w0 = 6.2831853071795862 * fcut / fsamp;
	c = cos(w0);
	for (k = n / 2; k >= 1; --k) {
		q = -0.5 / cos(3.1415926535897932 * (double) (2 * k + n - 1) / (double) (2 * n));
		r = sin(w0) / (2. * q);
		if (lo) {
			a1 = (1. - c) / (1. + r);
			a0 = +0.5 * a1;
		} else {
			a1 = -(1. + c) / (1. + r);
			a0 = -0.5 * a1;
		}
		a2 = a0;
		b1 = -2. * c / (1. + r);
		b2 = (1. - r) / (1. + r);
/* Process the buffer. */
		x = 0.;
		x1 = 0.;
		x2 = 0.;
		y = 0.;
		y1 = 0.;
		y2 = 0.;
		for (j = 0; j < l; ++j) {
			x2 = x1;
			x1 = x;
			x = s[j];
			y2 = y1;
			y1 = y;
			y = a0 * x + a1 * x1 + a2 * x2 - b1 * y1 - b2 * y2;
			s[j] = y;
		}
	}
	return 0;
} /* end of bfilt */