BilinearDesignLH.cpp

00001 //------------------------------------------------------------------------
00002 //  Design of Butterworth LPF and HPF using bilinear transform
00003 //
00004 //   2017/03/25, Copyright (c) 2017 MIKAMI, Naoki
00005 //------------------------------------------------------------------------
00006 
00007 #include "BilinearDesignLH.hpp"
00008 
00009 namespace Mikami
00010 {
00011     // Execute design
00012     //      input
00013     //          fc: Cutoff frequency
00014     //          pb: Passband (LPF or HPF)
00015     //      output
00016     //          c : Coefficients for cascade structure
00017     //          g : Gain factor for cascade structure
00018     void BilinearDesign::Execute(float fc, Type pb, Coefs c[], float& g)
00019     {
00020         Butterworth();
00021         Bilinear(fc);
00022         ToCascade(pb);
00023         GetGain(pb);
00024         GetCoefs(c, g);
00025     }
00026 
00027     // Get poles for Butterworth characteristics
00028     void BilinearDesign::Butterworth()
00029     {
00030         float pi_2order = PI_/(2.0f*ORDER_);
00031         for (int j=0; j<ORDER_/2; j++)  // Pole with imaginary part >= 0
00032         {
00033             float theta = (2.0f*j + 1.0f)*pi_2order;
00034             sP_[j] = Complex(-cosf(theta), sinf(theta));
00035         }
00036     }
00037 
00038     // Bilinear transform
00039     //      fc: Cutoff frequency
00040     void BilinearDesign::Bilinear(float fc)
00041     {
00042         float wc = tanf(fc*PI_FS_);
00043         for (int k=0; k<ORDER_/2; k++)
00044             zP_[k] = (1.0f + wc*sP_[k])/(1.0f - wc*sP_[k]);
00045     }
00046 
00047     // Convert to coefficients for cascade structure
00048     void BilinearDesign::ToCascade(Type pb)
00049     {
00050         for (int j=0; j<ORDER_/2; j++)
00051         {
00052             ck_[j].a1 = 2.0f*real(zP_[j]);          // a1m
00053             ck_[j].a2 = -norm(zP_[j]);              // a2m
00054             ck_[j].b1 = (pb == LPF) ? 2.0f : -2.0f; // b1m
00055             ck_[j].b2 = 1.0f;                       // b2m
00056         }
00057     }
00058 
00059     // Calculate gain factor
00060     void BilinearDesign::GetGain(Type pb)
00061     {
00062         float u = (pb == LPF) ? 1.0f : -1.0f;   // u: inverse of z
00063         gain_ = 1.0f;
00064         for (int k=0; k<ORDER_/2; k++)
00065             gain_ = gain_*(1.0f - (ck_[k].a1 + ck_[k].a2*u)*u)/
00066                           (1.0f + (ck_[k].b1 + ck_[k].b2*u)*u);
00067     }
00068 
00069     // Get coefficients
00070     void BilinearDesign::GetCoefs(Coefs c[], float& gain)
00071     {
00072         for (int k=0; k<ORDER_/2; k++) c[k] = ck_[k];
00073         gain = gain_;
00074     }
00075 }