File content as of revision 16:d60e5187fd31:

//------------------------------------------------------------------------
//  Design of Butterworth LPF and HPF using bilinear transform
//
//   2017/03/25, Copyright (c) 2017 MIKAMI, Naoki
//------------------------------------------------------------------------

#include "BilinearDesignLH.hpp"

namespace Mikami
{
    // Execute design
    //      input
    //          fc: Cutoff frequency
    //          pb: Passband (LPF or HPF)
    //      output
    //          c : Coefficients for cascade structure
    //          g : Gain factor for cascade structure
    void BilinearDesign::Execute(float fc, Type pb, Coefs c[], float& g)
    {
        Butterworth();
        Bilinear(fc);
        ToCascade(pb);
        GetGain(pb);
        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(Type pb)
    {
        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
            ck_[j].b2 = 1.0f;                       // b2m
        }
    }

    // Calculate gain factor
    void BilinearDesign::GetGain(Type pb)
    {
        float u = (pb == LPF) ? 1.0f : -1.0f;   // u: inverse of z
        gain_ = 1.0f;
        for (int k=0; k<ORDER_/2; k++)
            gain_ = gain_*(1.0f - (ck_[k].a1 + ck_[k].a2*u)*u)/
                          (1.0f + (ck_[k].b1 + ck_[k].b2*u)*u);
    }

    // Get coefficients
    void BilinearDesign::GetCoefs(Coefs c[], float& gain)
    {
        for (int k=0; k<ORDER_/2; k++) c[k] = ck_[k];
        gain = gain_;
    }
}