File content as of revision 0:33908268d9ea:

//------------------------------------------------------------------------------
//  Design of Butterworth LPF and HPF using bilinear transform
//
//   2014/06/29, Copyright (c) 2014 MIKAMI, Naoki
//------------------------------------------------------------------------------

#include "BilinearDesignLH.hpp"

namespace Mikami
{
    // Execute design
    //      input
    //          fc: Cutoff frequency
    //      output
    //          c : Coefficients for cascade structure
    //          g : Gain factor for cascade structure
    void BilinearDesign::Execute(float fc, Coefs c[], float& g)
    {
        Butterworth();
        Bilinear(fc);
        ToCascade();
        GetGain();
        GetCoefs(c, g);
    }

    // Get poles for Butterworth characteristics
    void BilinearDesign::Butterworth()
    {
        float pi_2order = PI_/(2.0f*ORDER_);
        for (int j=0; j<ORDER_/2; j++)  // Pole with imaginary part >= 0
        {
            float theta = (2.0f*j + 1.0f)*pi_2order;
            sP_[j] = Complex(-cosf(theta), sinf(theta));
        }
    }
    
    // Bilinear transform
    //      fc: Cutoff frequency
    void BilinearDesign::Bilinear(float fc)
    {
        float wc = tanf(fc*PI_FS_);
        for (int k=0; k<ORDER_/2; k++)
            zP_[k] = (1.0f + wc*sP_[k])/(1.0f - wc*sP_[k]);
    }
    
    // Convert to coefficients for cascade structure
    void BilinearDesign::ToCascade()
    {
        for (int j=0; j<ORDER_/2; j++)
        {
            ck_[j].a1 = 2.0f*real(zP_[j]);              // a1m
            ck_[j].a2 = -norm(zP_[j]);                  // a2m
            ck_[j].b1 = (PB_ == LPF) ? 2.0f : -2.0f;    // b1m
        }
    }

    // Calculate gain factor
    void BilinearDesign::GetGain(){
        float u = (PB_ == LPF) ? 1.0f : -1.0f;
        float g0 = 1.0f;
        for (int k=0; k<ORDER_/2; k++)
            g0 = g0*(1.0f - (ck_[k].a1 + ck_[k].a2*u)*u)/
                    (1.0f + (ck_[k].b1 + u)*u);
        gain_ = g0;
    }
    
    // Get coefficients
    void BilinearDesign::GetCoefs(Coefs c[], float& gain)
    {
        for (int k=0; k<ORDER_/2; k++)
        {
            c[k].a1 = ck_[k].a1;
            c[k].a2 = ck_[k].a2;
            c[k].b1 = ck_[k].b1;
        }
        gain = gain_;
    }
}