bilinear.cpp

来自「Digital filter designer s handbook C++ c」· C++ 代码 · 共 167 行

//
//  File = bilinear.cpp
//                     

 #include <stdlib.h> 
 #include <iostream.h>
 #include <fstream.h> 
 #include "misdefs.h"
 #include "typedefs.h" 
 #include "d_cmplx.h"
 #include "bilinear.h"
 extern ofstream DebugFile;
 
//====================================================
//
//----------------------------------------------------
IirFilterDesign* BilinearTransf( 
                         FilterTransFunc* analog_filter,
                         double sampling_interval)

{
 int max_poles;
 int num_poles, num_zeros;
 int j, m, n;
 int max_coeff, num_numer_coeff;
 double h_const, h_sub_zero, denom_mu_zero;
 double_complex *pole, *zero;
 double_complex *mu;
 double_complex alpha;
 double_complex beta;
 double_complex gamma;
 double_complex delta;
 double_complex eta;
 double_complex work;
 double_complex c_two;
 double *a, *b;
 IirFilterDesign* iir_filter;
 
 pole = analog_filter->GetPoles(&num_poles);
 zero = analog_filter->GetZeros(&num_zeros);
 h_sub_zero = analog_filter->GetHSubZero();
 DebugFile << "num analog poles = " << num_poles << endl;
 DebugFile << "num analog zeros = " << num_zeros << endl;

 if(num_poles > num_zeros)
   { max_poles = num_poles; }
 else
   { max_poles = num_zeros; }
   
 //--------------------------------------------
 // allocate and initialize working storage
 
 mu = new double_complex[max_poles+1];
 a = new double[max_poles+1];
 b = new double[max_poles+1];
 
 for(j=0; j<=max_poles; j++)
   {
    mu[j] = double_complex(0.0, 0.0);
    a[j] = 0.0;
    b[j] = 0.0;
   }           
   
 //-------------------------------------------
 // compute constant gain factor
 
 h_const = 1.0;
 work = double_complex(1.0, 0.0);
 c_two = double_complex(2.0, 0.0);
 
 for(n=1; n<=num_poles; n++)
   {
    work = work * (c_two - (sampling_interval * pole[n]));
    h_const = h_const * sampling_interval;
   } 
 DebugFile << "T**2 = " << h_const << endl;                         
 h_const = h_sub_zero * h_const / real(work);
 
 DebugFile << "in BilinearTransf, h_const = "
           << h_const << endl;
 DebugFile << "work = " << work << endl;
 
 //--------------------------------------------------
 // compute denominator coefficients
 
 mu[0] = double_complex(1.0, 0.0);
 
 for(n=1; n<=num_poles; n++)
   {
    DebugFile << "in BilinearTransf, pole [" << n
              << "] = " << pole[n] << endl;
    gamma = double_complex( (2.0/sampling_interval), 0.0) - pole[n];
    delta = double_complex( (-2.0/sampling_interval), 0.0) - pole[n];
    DebugFile << "gamma = " << gamma << endl;
    DebugFile << "delta = " << delta << endl;
    
    for(j=n; j>=1; j--)
      {
       mu[j] = gamma * mu[j] + (delta * mu[j-1]);
      }
    mu[0] = gamma * mu[0];
   } 
 DebugFile << "for denom, mu[0] = " << mu[0] << endl;  
 denom_mu_zero = real(mu[0]);                                
 for( j=1; j<=num_poles; j++)
   {
    //a[j] = -h_const * real(mu[j]) /h_sub_zero;
    a[j] = -1.0 * real(mu[j])/denom_mu_zero;
    //a[j] = -1.0 * real(mu[num_poles+1 - j]);
    DebugFile << "a[" << j << "] = " << a[j] << endl;
    DebugFile << "imag(mu[" << j << "]) = "
              << imag(mu[j]) << endl;
   }
 //-----------------------------------------------------
 //  compute numerator coeffcients
 
 mu[0] = double_complex(1.0, 0.0);
 for(n=1; n<=max_poles; n++)
   {
    mu[n] = double_complex(0.0, 0.0);
   }                                 
   
max_coeff = 0;
 
 //- - - - - - - - - - - - - - - - - - - - -
 //  compute (1+z**(-1)) ** (N-M)
 
 for(m=1; m<=(num_poles-num_zeros); m++)
   {
    max_coeff++;
    for(j=max_coeff; j>=1; j--)
      {
       mu[j] = mu[j] + mu[j-1];
      }
   }
 for(m=1; m<=num_zeros; m++)
   {
    max_coeff++;
    DebugFile << "zero[" << m << "] = " << zero[m] << endl;
    alpha = double_complex( (2.0/sampling_interval), 0.0) - zero[m];
    //beta = double_complex( (-2.0/sampling_interval), 0.0) - zero[m]; 
    beta = double_complex( (-2.0/sampling_interval), 0.0) - zero[m]; 
    
    for(j=max_coeff; j>=1; j--)
      {
       mu[j] = alpha * mu[j] + (beta * mu[j-1]);
      }
    mu[0] = alpha * mu[0];
   }
 num_numer_coeff = max_coeff+1;
 for(j=0; j<num_numer_coeff; j++)
   {
    //b[j] = h_const * real(mu[j]);
    b[j] = h_sub_zero * real(mu[j])/denom_mu_zero;
    DebugFile << "b[" << j << "] = " << b[j] << endl;
    DebugFile << "imag(mu[" << j << "]) = " 
              << imag(mu[j]) << endl;
   }                             
 
 delete []mu; 
 iir_filter = new IirFilterDesign( num_numer_coeff, 
                                   num_poles,
                                   b, a);
 iir_filter->SetSamplingInterval(sampling_interval);                 
 return(iir_filter);
}