BilinearDesignLH.cpp

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