?

Committer:
phungductung
Date:
Tue Jun 04 21:58:08 2019 +0000
Revision:
0:cacfc9e25452
?;

Who changed what in which revision?

UserRevisionLine numberNew contents of line
phungductung 0:cacfc9e25452 1 //------------------------------------------------------------------------------
phungductung 0:cacfc9e25452 2 // Design of Butterworth LPF and HPF using bilinear transform
phungductung 0:cacfc9e25452 3 //
phungductung 0:cacfc9e25452 4 // 2016/03/31, Copyright (c) 2016 MIKAMI, Naoki
phungductung 0:cacfc9e25452 5 //------------------------------------------------------------------------------
phungductung 0:cacfc9e25452 6
phungductung 0:cacfc9e25452 7 #include "BilinearDesignLH.hpp"
phungductung 0:cacfc9e25452 8
phungductung 0:cacfc9e25452 9 namespace Mikami
phungductung 0:cacfc9e25452 10 {
phungductung 0:cacfc9e25452 11 // Execute design
phungductung 0:cacfc9e25452 12 // input
phungductung 0:cacfc9e25452 13 // fc: Cutoff frequency
phungductung 0:cacfc9e25452 14 // pb: Passband (LPF or HPF)
phungductung 0:cacfc9e25452 15 // output
phungductung 0:cacfc9e25452 16 // c : Coefficients for cascade structure
phungductung 0:cacfc9e25452 17 // g : Gain factor for cascade structure
phungductung 0:cacfc9e25452 18 void BilinearDesign::Execute(float fc, Type pb, Coefs c[], float& g)
phungductung 0:cacfc9e25452 19 {
phungductung 0:cacfc9e25452 20 Butterworth();
phungductung 0:cacfc9e25452 21 Bilinear(fc);
phungductung 0:cacfc9e25452 22 ToCascade(pb);
phungductung 0:cacfc9e25452 23 GetGain(pb);
phungductung 0:cacfc9e25452 24 GetCoefs(c, g);
phungductung 0:cacfc9e25452 25 }
phungductung 0:cacfc9e25452 26
phungductung 0:cacfc9e25452 27 // Get poles for Butterworth characteristics
phungductung 0:cacfc9e25452 28 void BilinearDesign::Butterworth()
phungductung 0:cacfc9e25452 29 {
phungductung 0:cacfc9e25452 30 float pi_2order = PI_/(2.0f*ORDER_);
phungductung 0:cacfc9e25452 31 for (int j=0; j<ORDER_/2; j++) // Pole with imaginary part >= 0
phungductung 0:cacfc9e25452 32 {
phungductung 0:cacfc9e25452 33 float theta = (2.0f*j + 1.0f)*pi_2order;
phungductung 0:cacfc9e25452 34 sP_[j] = Complex(-cosf(theta), sinf(theta));
phungductung 0:cacfc9e25452 35 }
phungductung 0:cacfc9e25452 36 }
phungductung 0:cacfc9e25452 37
phungductung 0:cacfc9e25452 38 // Bilinear transform
phungductung 0:cacfc9e25452 39 // fc: Cutoff frequency
phungductung 0:cacfc9e25452 40 void BilinearDesign::Bilinear(float fc)
phungductung 0:cacfc9e25452 41 {
phungductung 0:cacfc9e25452 42 float wc = tanf(fc*PI_FS_);
phungductung 0:cacfc9e25452 43 for (int k=0; k<ORDER_/2; k++)
phungductung 0:cacfc9e25452 44 zP_[k] = (1.0f + wc*sP_[k])/(1.0f - wc*sP_[k]);
phungductung 0:cacfc9e25452 45 }
phungductung 0:cacfc9e25452 46
phungductung 0:cacfc9e25452 47 // Convert to coefficients for cascade structure
phungductung 0:cacfc9e25452 48 void BilinearDesign::ToCascade(Type pb)
phungductung 0:cacfc9e25452 49 {
phungductung 0:cacfc9e25452 50 for (int j=0; j<ORDER_/2; j++)
phungductung 0:cacfc9e25452 51 {
phungductung 0:cacfc9e25452 52 ck_[j].a1 = 2.0f*real(zP_[j]); // a1m
phungductung 0:cacfc9e25452 53 ck_[j].a2 = -norm(zP_[j]); // a2m
phungductung 0:cacfc9e25452 54 ck_[j].b1 = (pb == LPF) ? 2.0f : -2.0f; // b1m
phungductung 0:cacfc9e25452 55 ck_[j].b2 = 1.0f; // b2m
phungductung 0:cacfc9e25452 56 }
phungductung 0:cacfc9e25452 57 }
phungductung 0:cacfc9e25452 58
phungductung 0:cacfc9e25452 59 // Calculate gain factor
phungductung 0:cacfc9e25452 60 void BilinearDesign::GetGain(Type pb)
phungductung 0:cacfc9e25452 61 {
phungductung 0:cacfc9e25452 62 float u = (pb == LPF) ? 1.0f : -1.0f;
phungductung 0:cacfc9e25452 63 float g0 = 1.0f;
phungductung 0:cacfc9e25452 64 for (int k=0; k<ORDER_/2; k++)
phungductung 0:cacfc9e25452 65 g0 = g0*(1.0f - (ck_[k].a1 + ck_[k].a2*u)*u)/
phungductung 0:cacfc9e25452 66 (1.0f + (ck_[k].b1 + ck_[k].b2*u)*u);
phungductung 0:cacfc9e25452 67 gain_ = g0;
phungductung 0:cacfc9e25452 68 }
phungductung 0:cacfc9e25452 69
phungductung 0:cacfc9e25452 70 // Get coefficients
phungductung 0:cacfc9e25452 71 void BilinearDesign::GetCoefs(Coefs c[], float& gain)
phungductung 0:cacfc9e25452 72 {
phungductung 0:cacfc9e25452 73 for (int k=0; k<ORDER_/2; k++) c[k] = ck_[k];
phungductung 0:cacfc9e25452 74 gain = gain_;
phungductung 0:cacfc9e25452 75 }
phungductung 0:cacfc9e25452 76 }