anlg_filt_intg.cpp
上传用户:jtjnyq9001
上传日期:2014-11-21
资源大小:3974k
文件大小:9k
源码类别:

3G开发

开发平台:

Visual C++

  1. //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2. //
  3. //  File = anlg_filt_intg.cpp
  4. //
  5. //  The class AnalogFilterByInteg serves as the base class
  6. //  for all the "traditional" analog filter types
  7. //  e.g. Butterworth, Chebyshev, etc. that are modeled 
  8. //  using numerical integration
  9. //
  10. #include <math.h>
  11. #include "misdefs.h"
  12. #include "parmfile.h"
  13. #include "model_graph.h"
  14. //#include "iir_filt.h"
  15. #include "anlg_filt_intg.h"
  16. #include "bilin_transf.h"
  17. #include "iir_comp_resp.h"
  18. extern ParmFile *ParmInput;
  19. #ifdef _DEBUG
  20.   extern ofstream *DebugFile;
  21. #endif
  22. //======================================================
  23. //  constructor
  24. AnalogFilterByInteg::AnalogFilterByInteg( char *instance_name,
  25.                         PracSimModel *outer_model,
  26.                         Signal<float> *in_sig,
  27.                         Signal<float> *out_sig )
  28.            :PracSimModel( instance_name,
  29.                           outer_model )
  30. {
  31.   OPEN_PARM_BLOCK;
  32.   In_Sig = in_sig;
  33.   Out_Sig = out_sig;
  34.   Denorm_Proto_Filt = NULL;
  35.   Lowpass_Proto_Filt = NULL;
  36.   A_Coefs = NULL;
  37.   B_Coefs = NULL;
  38.   //In_Mem = NULL;
  39.   //Out_Mem = NULL;
  40.   Filt_Order = 0;
  41.   Estim_Delay = 0;
  42.   GET_BOOL_PARM(Bypass_Enabled);
  43.   if(Bypass_Enabled)
  44.     {
  45.     Filt_Band_Config = LOWPASS_FILT_BAND_CONFIG;
  46.     Filt_Order = 0;
  47.     Norm_Hz_Pass_Edge = 0.0;
  48.     Norm_Hz_Pass_Edge_2 = 0.0;
  49.     Lowpass_Proto_Filt = NULL;
  50.     Denorm_Proto_Filt = NULL;
  51.     }
  52.   else
  53.     {
  54.     GET_BOOL_PARM(Using_Biquads);
  55.     GET_INT_PARM(Filt_Order);
  56.     Filt_Band_Config = GetFiltBandConfigParm("Filt_Band_Config");
  57.     GET_DOUBLE_PARM(Integ_Alpha);
  58.     switch (Filt_Band_Config)
  59.       {
  60.       case LOWPASS_FILT_BAND_CONFIG:
  61.         GET_DOUBLE_PARM(Norm_Hz_Pass_Edge);
  62.         Norm_Hz_Pass_Edge_2 = 0;
  63.         Prototype_Order = Filt_Order;
  64.         break;
  65.       case HIGHPASS_FILT_BAND_CONFIG:
  66.         GET_DOUBLE_PARM(Norm_Hz_Pass_Edge);
  67.         Norm_Hz_Pass_Edge_2 = 0;
  68.         Prototype_Order = Filt_Order;
  69.         break;
  70.       case BANDPASS_FILT_BAND_CONFIG:
  71.       case BANDSTOP_FILT_BAND_CONFIG:
  72.         {
  73.         GET_DOUBLE_PARM(Norm_Hz_Pass_Edge);
  74.         GET_DOUBLE_PARM(Norm_Hz_Pass_Edge_2);
  75.         if( Norm_Hz_Pass_Edge >= Norm_Hz_Pass_Edge_2)
  76.           {
  77.           // error
  78.           }
  79.         if( Filt_Order % 2 != 0 )
  80.           {
  81.           // another error
  82.           }
  83.         Prototype_Order = Filt_Order/2;
  84.         break;
  85.         } // end of cases for bandpass and bandstop
  86.       } // end of switch on Filter_Band_Config
  87.     GET_BOOL_PARM(Resp_Plot_Enabled);
  88.     if(Resp_Plot_Enabled)
  89.       {
  90.       GET_BOOL_PARM(Db_Scale_Enabled);
  91.       }
  92.     }
  93.   MAKE_INPUT(In_Sig);
  94.   MAKE_OUTPUT(Out_Sig);
  95.   Cumul_Samp_Count = 0;
  96.   return;
  97. }
  98. //=============================================================================
  99. // destructor
  100. AnalogFilterByInteg::~AnalogFilterByInteg(){};
  101. //=============================================================================
  102. void AnalogFilterByInteg::Initialize(int proc_block_size, double samp_intvl)
  103. {
  104.   Proc_Block_Size = proc_block_size;
  105.   Samp_Intvl = samp_intvl;
  106.   Init_Kernel();
  107. }
  108. //=============================================================================
  109. void AnalogFilterByInteg::Initialize()
  110. {
  111.   Samp_Intvl = In_Sig->GetSampIntvl();
  112.   Proc_Block_Size = In_Sig->GetBlockSize();
  113.   Init_Kernel();
  114. }
  115. //=============================================================================
  116. void AnalogFilterByInteg::Init_Kernel()
  117. {
  118.   Warped_Lower_Passband_Edge = Norm_Hz_Pass_Edge;
  119.   Warped_Upper_Passband_Edge = Norm_Hz_Pass_Edge_2;
  120.   //Warped_Lower_Passband_Edge =
  121.   //          tan( PI * Samp_Intvl * Norm_Hz_Pass_Edge )/
  122.   //          (PI * Samp_Intvl);
  123.   //Warped_Upper_Passband_Edge =
  124.   //          tan( PI * Samp_Intvl * Norm_Hz_Pass_Edge_2 )/
  125.   //          (PI * Samp_Intvl);
  126.   Denorm_Proto_Filt = new DenormalizedPrototype(  Lowpass_Proto_Filt,
  127.                                                   Filt_Band_Config,
  128.                                                   Warped_Lower_Passband_Edge,
  129.                                                   Warped_Upper_Passband_Edge);
  130.    if(Using_Biquads)
  131.    {
  132.       Integrator_0 = new NumericInteg*[Num_Biquads];
  133.       Integrator_1 = new NumericInteg*[Num_Biquads];
  134.       W0_Prime = new double[Num_Biquads];
  135.       W1_Prime = new double[Num_Biquads];
  136.       W2_Prime = new double[Num_Biquads];
  137.       B0_Coef = new double[Num_Biquads];
  138.       B1_Coef = new double[Num_Biquads];
  139.       A0_Coef = new double[Num_Biquads];
  140.       A1_Coef = new double[Num_Biquads];
  141.       A2_Coef = new double[Num_Biquads];
  142.    }
  143.    else
  144.    {  
  145.       Integrator = new NumericInteg*[Filt_Order];
  146.       Y_Prime = new double[Filt_Order+1];
  147.       B_Coefs = new double[Filt_Order+1];
  148.       A_Coefs = new double[Filt_Order];
  149.    }
  150.   H_Zero = Denorm_Proto_Filt->GetHZero();
  151.   //H_Zero = Lowpass_Proto_Filt->GetHZero();
  152.   #ifdef _DEBUG
  153.     *DebugFile << "H_Zero = " << H_Zero << endl;
  154.   #endif
  155.   int k;
  156.   *DebugFile << "in AnalofFilterByInteg, Lowpass B_coeffs:" << endl;
  157.   Lowpass_Proto_Filt->GetDenomPoly(false);
  158.   *DebugFile << "check point A" << endl;
  159.   for( k=0; k<Filt_Order; k++ )
  160.     {
  161.     *DebugFile << "k = " << k << endl;
  162.     Integrator[k] = new NumericInteg(Samp_Intvl, Integ_Alpha);
  163.     B_Coefs[k] = -(Denorm_Proto_Filt->GetDenomPoly(true).GetCoefficient(k));
  164.     //B_Coefs[k] = -(Lowpass_Proto_Filt->GetDenomPoly(false).GetCoefficient(k));
  165.     Y_Prime[k] = 0.0;
  166.     #ifdef _DEBUG
  167.       *DebugFile << "in AnalogFilterByInteg, B_Coefs["
  168.                 << k << "] = " << B_Coefs[k] << endl;
  169.     #endif
  170.     }
  171.   
  172.   #ifdef _DEBUG
  173.     int idbg;
  174. //    for(idbg=1; idbg<=Filt_Order; idbg++)
  175. //      {
  176. //      *DebugFile << "A_Coefs[" << idbg << "] = " << A_Coefs[idbg] << endl;
  177. //      }
  178.     for(idbg=0; idbg<=Filt_Order; idbg++)
  179.       {
  180.       *DebugFile << "B_Coefs[" << idbg << "] = " << B_Coefs[idbg] << endl;
  181.       }
  182.   #endif
  183.    Num_Zeros = Denorm_Proto_Filt->GetNumZeros();
  184.    if(Num_Zeros == 0)
  185.    {
  186.       for( k=1; k<Filt_Order; k++)
  187.       {
  188.          A_Coefs[k] = 0.0;
  189.       }
  190.       A_Coefs[0] = 1.0;
  191.    }
  192.    else
  193.    {
  194.       for( k=0; k<=Num_Zeros; k++)
  195.       {
  196.          A_Coefs[k] = (Denorm_Proto_Filt->GetNumerPoly(true).GetCoefficient(k));
  197.       }
  198.    }
  199. }
  200. //===========================================
  201. //  Method to execute the model
  202. int AnalogFilterByInteg::Execute()
  203. {
  204.    float *in_sig_ptr;
  205.    float *out_sig_ptr;
  206.    double sum, output;
  207.    int samp_idx;
  208.    int sec, k;
  209.    int proc_block_size;
  210.    in_sig_ptr = GET_INPUT_PTR(In_Sig);
  211.    out_sig_ptr = GET_OUTPUT_PTR(Out_Sig);
  212.   proc_block_size = In_Sig->GetValidBlockSize();
  213.   Out_Sig->SetValidBlockSize(proc_block_size);
  214.    if(Bypass_Enabled)
  215.    { // copy input signal to output signal
  216.       for(samp_idx=0; samp_idx<Proc_Block_Size; samp_idx++)
  217.       {
  218.          *out_sig_ptr = *in_sig_ptr;
  219.          out_sig_ptr++;
  220.          in_sig_ptr++;
  221.       }
  222.    }
  223.    else if(Using_Biquads)
  224.    {
  225.       for(samp_idx=0; samp_idx<Proc_Block_Size; samp_idx++)
  226.       {
  227.          sum = (*in_sig_ptr) * H_Zero;
  228.          for(sec=0; sec<Num_Biquads; sec++)
  229.          {
  230.             sum += W0_Prime[sec] * B0_Coef[sec];
  231.             sum += W1_Prime[sec] * B1_Coef[sec];
  232.             W2_Prime[sec] = sum;
  233.             W1_Prime[sec] = Integrator_1[sec]->Integrate(W2_Prime[sec]);
  234.             W0_Prime[sec] = Integrator_0[sec]->Integrate(W1_Prime[sec]);
  235.          }
  236.          sum = A0_Coef[sec] * W0_Prime[sec]
  237.                   + A1_Coef[sec] * W1_Prime[sec]
  238.                   + A2_Coef[sec] * W2_Prime[sec];
  239.          *out_sig_ptr = sum;
  240.          out_sig_ptr++;
  241.          in_sig_ptr++;
  242.          Cumul_Samp_Count++;
  243.       }
  244.    }
  245.    else // using individual coefficients
  246.    {
  247.       proc_block_size = Proc_Block_Size;
  248.       //  perform filter summations
  249.       for(samp_idx=0; samp_idx<proc_block_size; samp_idx++)
  250.       {
  251.          sum = (*in_sig_ptr) * H_Zero;
  252.          for(k=0; k<Filt_Order; k++)
  253.          {
  254.             sum += (Y_Prime[k] * B_Coefs[k]);
  255.          }
  256.          Y_Prime[Filt_Order] = sum;
  257.          output = 0.0;
  258.          for( k=Filt_Order-1; k>=0; k--)
  259.          {
  260.             Y_Prime[k] = ((Integrator[k])->Integrate(Y_Prime[k+1]));
  261.             output += (Y_Prime[k] * A_Coefs[k]);
  262.          }
  263.          //*out_sig_ptr = H_Zero * output;
  264.          *out_sig_ptr = output;
  265.          out_sig_ptr++;
  266.          in_sig_ptr++;
  267.          Cumul_Samp_Count++;
  268.       } // end of loop over samp_idx
  269.    } // end of else clause on if(Bypass_Enabled) control structure
  270.   return(_MES_AOK);
  271. }