anlg_sim.cpp

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

//
//  file = anlg_sim.cpp
//

 #include <stdlib.h> 
 #include <iostream.h> 
 #include <fstream.h>
 #include <math.h> 
 #include "polefilt.h"
 #include "pzfilt.h"
 #include "anlgfilt.h"
 #include "cmpxpoly.h"
 #include "poly.h"
 #include "filtfunc.h"
 #include "buttfunc.h"
 #include "chebfunc.h"
 #include "elipfunc.h"
  
 ofstream DebugFile("analog.out", ios::out);

 main()
 {
  double output_sample;
  double delta_t;
  float total_time;
  double passband_ripple, stopband_ripple;
  double passband_edge, stopband_edge;
  long int k;
  int order, plotter_decimation; 
  int type_of_filter; 
  int type_of_test_signal;
  int ripple_bw_norm;
  AnalogFilter* filter;
  
  CmplxPolynomial cmplx_denom_poly;
  Polynomial denom_poly;
  CmplxPolynomial cmplx_numer_poly;
  Polynomial numer_poly;
  FilterTransFunc* filter_function;
  ofstream* OutputFile;
    
  cout << "This program simulates the operation\n"
       << "of a 'classical' analog filter\n" << endl;
  cout << "Desired type of filter:\n"
       << "  1 = Butterworth\n"
       << "  2 = Chebyshev\n"
       << "  3 = Chebyshev Type II (N/A)\n"
       << "  4 = Elliptical\n"
       << "  5 = Bessel(N/A)\n"
       << endl;
  cin >> type_of_filter;
  
  cout << "Desired order of filter: " << endl;
  cin >> order;
  
  if((type_of_filter == 2) || (type_of_filter == 4))
    {
     cout << "Desired passband ripple: " << endl;
     cin >> passband_ripple;
    }
  
  if((type_of_filter == 3) || (type_of_filter == 4))
    {
     cout << "Desired stopband ripple: " << endl;
     cin >> stopband_ripple;
    }
  
  cout << "Desired passband edge: " << endl;
  cin >> passband_edge;
  
  if( type_of_filter == 2 )
    {
     cout << "Desired type of normalization?\n"
          << "  0 = 3 dB bandwidth\n"
          << "  1 = ripple bandwidth" << endl;
     cin >> ripple_bw_norm;
    }
  
  if(type_of_filter == 4)
    {
     cout << "Desired stopband edge: " << endl;
     cin >> stopband_edge;
    }
  
  cout << "Desired type of test signal:\n"
       << "   1 = impulse response\n"
       << "   2 = step response" << endl;
  cin >> type_of_test_signal;
  
  try_again:
  cout << "Desired time increment: " << endl;
  cin >> delta_t;
  
  if(type_of_filter ==4)
    {
    if( (1./delta_t) < (2.0*stopband_edge) )
      {
      cout << "for this filter, time increment must be < "
           << 0.5/stopband_edge << endl;
      goto try_again;
      }
    } 
  else
    {
    if( (1./delta_t) < (5.0*passband_edge) )
      {
      cout << "for this filter, time increment must be < "
           << 0.2/passband_edge << endl;
      goto try_again;
      }
    }  
    
  cout << "Downsampling rate for plotter output" << endl;
  cin >> plotter_decimation;
  
  cout << "Desired duration of simulation: " << endl;
  cin >> total_time;

  switch (type_of_filter)
    {
    case 1:          // Butterworth
      filter_function = new ButterworthTransFunc(order);
      filter_function->LowpassDenorm(passband_edge);
      break;
    case 2:
      filter_function = new ChebyshevTransFunc( order, 
                                               (float)passband_ripple,
                                               ripple_bw_norm);
      filter_function->LowpassDenorm(passband_edge);
      break;
    case 3:
      cout << "Chebyshev type II not supported yet" << endl;
      return 0;
    case 4:
      int upper_summation_limit = 5;
      filter_function = new EllipticalTransFunc( order, 
												                         passband_ripple,
                                                 stopband_ripple,
                                                 passband_edge,
                                                 stopband_edge,
                                                 upper_summation_limit);
      filter_function->DumpBiquads(&DebugFile);
      filter_function->LowpassDenorm(
                       sqrt(passband_edge * stopband_edge));
      break;
    
    }

  if(type_of_filter == 4)
    {
     filter = new AnalogPoleZeroFilter( filter_function->GetNumerPoly(),
                                        filter_function->GetDenomPoly(),
                                        filter_function->GetHSubZero(),
                                        delta_t);
    }
  else
    {
     filter = new AnalogAllPoleFilt( filter_function->GetDenomPoly(),
                                     filter_function->GetHSubZero(),
                                     delta_t);
    }

 
  long int k_stop = (long int)(total_time/delta_t);
  cout << "k_stop = " << k_stop << endl;  

  switch (type_of_test_signal)
    {
    case 1:   // impulse response
      OutputFile = new ofstream("imp_resp.txt",ios::out);
      k=0; 
      output_sample = filter->Run( 1.0/delta_t );
      (*OutputFile) << k << ", " << output_sample << endl;
  
      for( k=1; k<k_stop; k++)
        {
         output_sample = filter->Run(0.0);
           if( (k % plotter_decimation)==0 )
             {
              (*OutputFile) << k*delta_t << ", " << output_sample << endl;
             }
        } 
      break;
    case 2:   //step response
      OutputFile = new ofstream("stp_resp.txt",ios::out);
      k=0; 
      output_sample = filter->Run(1.0);
      (*OutputFile) << k << ", " << output_sample << endl;
  
      for( k=1; k<k_stop; k++)
        {
         output_sample = filter->Run(1.0);
         if(k==20)
          {
          cout << "k = " << k << endl;
          cout << "plotter_decimation = " << plotter_decimation << endl;
          cout << "k%plotter_decimation = " << k%plotter_decimation << endl;
          }
         if( (k % plotter_decimation) == 0 )
           {
            (*OutputFile) << k*delta_t << ", " << output_sample << endl;  
           }
        } 
      break;
                
    }
     
  OutputFile->close();
  return 0;
 }