Butterworth bandpass filter biquad coefficients

butter.cxx

#define _USE_MATH_DEFINES

#include <cmath>
#include <complex>
#include <iostream>
#include <vector>

using namespace std;

// Based on the excellent article from Neil Robertson 
// https://www.dsprelated.com/showarticle/1257.php 

int main() {
  int N = 3; // Order of prototype LPF (number of biquads in BPF)
  double flo = 20-2;
  double fup = 20+2;
  double fs = 100; // Sampling rate
  
  std::vector<std::complex<double>> pa_lp(N); // Poles of prototype LPF (s-domain)
  std::vector<std::complex<double>> pa_bp(N); // Poles of BPF (s-domain)
  std::vector<std::complex<double>> pz_bp(N); // Poles of BPF (z-domain)
  
  // Poles of the prototype low-pass Butterworth filter (s-domain)
  for (int k = 0; k < N; k++) {
    double theta= ((2.0*k + 1.0)*M_PI)/(2.0*N);
    cout << "Theta: " << theta << endl;
    pa_lp[k].real(-sin(theta));
    pa_lp[k].imag(cos(theta));
    cout << pa_lp[k] << endl;
  }
  
  // Pre-warp frequencies to compensate for the non-linearity of the bilinear transform
  double Flo = fs/M_PI * tan(M_PI*flo/fs);
  double Fup = fs/M_PI * tan(M_PI*fup/fs);
  double BW = Fup-Flo; // Hz -3 dB bandwidth
  double F0 = sqrt(Flo*Fup); // Hz geometric mean frequency
  
  // Transform the prototype LPF poles to BPF poles, still in s-domain
  // Conjugates are not computed
  for (int k = 0; k < N; k++) {
    std::complex<double> alpha = (BW/(2.0 * F0)) * pa_lp[k];
    std::complex<double> beta = std::sqrt(std::complex<double>(1.0, 0.0) - std::pow((BW/(2.0 * F0)) * pa_lp[k], 2));
    beta *= std::complex<double>(0.0f, 1.0f);
    pa_bp[k] = (alpha + beta) * 2.0 * M_PI * F0;
  }

  // Bilinear transform poles from s-domain to z-domain
  for (int k = 0; k < N; k++) {
    pz_bp[k] = (pa_bp[k]/(2.0*fs) + 1.0) / (-pa_bp[k]/(2.0*fs) + 1.0);
  }

  // Denominator coeffs
  // Each biquad takes care of a complex conjugate pair of poles.
  // Coefficients a1, a2 are derived from: (z - (a + jb))(z - (a - jb))

  double f0 = std::sqrt(flo * fup);
  for (int k = 0; k < N; k++) {
    double a1 = -2.0*pz_bp[k].real();
    double a2 = std::norm(pz_bp[k]);
    std::cout << "a1: " << a1 << " a2: " << a2 << std::endl;

    // Biquad gain
    double fn = (2.0 * M_PI * f0) / fs; // Normalised frequency
    std::complex<double> z(cos(fn), sin(fn));
    std::complex<double> H = (std::pow(z, 2) - 1.0) / (std::pow(z, 2) + a1 * z + a2);
    cout << "k: " << 1.0 / abs(H) << endl;
  }
   
  return 0;  
}