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

3G开发

开发平台:

Visual C++

  1. //
  2. //  File = bart_pdgm_wind.cpp
  3. //
  4. #include <stdlib.h>
  5. #include "parmfile.h"
  6. #include "model_graph.h"
  7. #include "bart_pdgm_wind.h"
  8. #include "trianglr.h"
  9. #include "hamming.h"
  10. #include "hann.h"
  11. #include "fft_T.h"
  12. #include "dump_spect.h"
  13. #ifdef _DEBUG
  14.   extern ofstream *DebugFile;
  15. #endif
  16. #define _NO_ZERO_ENDS 0
  17. #define _ZERO_ENDS 1
  18. extern ParmFile *ParmInput;
  19. extern int PassNumber;
  20. //======================================================
  21. template <class T>
  22. BartlettPeriodogramWindowed<T>::
  23.             BartlettPeriodogramWindowed( 
  24.                               char* instance_name,
  25.                               PracSimModel* outer_model,
  26.                               Signal<T>* in_sig )
  27.                 :PracSimModel( instance_name,
  28.                                outer_model )
  29. {
  30.    int is;
  31.    MODEL_NAME(BartlettPeriodogram);
  32.    OPEN_PARM_BLOCK;
  33.    GET_INT_PARM(Seg_Len);
  34.    GET_INT_PARM(Fft_Len);
  35.    GET_INT_PARM(Hold_Off);
  36.    GET_INT_PARM(Num_Segs_To_Avg);
  37.    GET_DOUBLE_PARM(Freq_Norm_Factor);
  38.    GET_BOOL_PARM(Output_In_Decibels);
  39.    GET_BOOL_PARM(Plot_Two_Sided);
  40.    GET_BOOL_PARM(Halt_When_Completed);
  41.    GET_BOOL_PARM(Using_Window);
  42.    Psd_File_Name = new char[64];
  43.    strcpy(Psd_File_Name, "");
  44.    GET_STRING_PARM(Psd_File_Name);
  45.    if(Using_Window){
  46.       Window_Shape = 
  47.          GetWindowShapeParm("Window_Shape");
  48.       switch(Window_Shape){
  49.       case WINDOW_SHAPE_TRIANGULAR:
  50.          Data_Window = new TriangularWindow( Seg_Len,
  51.                                        _NO_ZERO_ENDS );
  52.          break;
  53.       case WINDOW_SHAPE_HAMMING:
  54.          Data_Window = new HammingWindow( Seg_Len );
  55.          break;
  56.       case WINDOW_SHAPE_HANN:
  57.          Data_Window = new HannWindow( Seg_Len, 
  58.                                        _NO_ZERO_ENDS );
  59.          break;
  60.       }
  61.       Window_Taps = Data_Window->GetDataWindow();
  62.       Window_Power = 0.0;
  63.       for(is=0; is<Seg_Len; is++){
  64.          Window_Power +=Window_Taps[is]*Window_Taps[is];
  65.       }
  66.       Window_Power /= double(Seg_Len);
  67.       double window_scale = sqrt(Window_Power);
  68.       for(is=0; is<Seg_Len; is++){
  69.          Window_Taps[is] /= window_scale;
  70.       }
  71.    }
  72.    In_Sig = in_sig;
  73.    MAKE_INPUT(In_Sig);
  74.    Time_Seg = new T[Seg_Len];
  75.    Win_Time_Seg = new T[Seg_Len];
  76.    Sample_Spectrum = new double[Fft_Len];
  77.    Freq_Seg = new std::complex<double>[Fft_Len];
  78.    for(is=0; is<Fft_Len; is++){
  79.       Sample_Spectrum[is] = 0.0;
  80.    }
  81.    Psd_File = new ofstream(Psd_File_Name, ios::out);
  82.    Processing_Completed = false;
  83. }
  84. //======================================================
  85. template <class T>
  86. BartlettPeriodogramWindowed<T>::~BartlettPeriodogramWindowed( void ){ };
  87. //======================================================
  88. template <class T>
  89. void BartlettPeriodogramWindowed<T>::Initialize(void)
  90. {
  91.    Segs_In_Est = 0;
  92.    Samps_Needed = Seg_Len;
  93.    Block_Size = In_Sig->GetBlockSize();
  94.    Samp_Intvl = In_Sig->GetSampIntvl();
  95.    Delta_F = 1.0/(Samp_Intvl*Fft_Len);
  96. };
  97. //======================================================
  98. template <class T>
  99. int BartlettPeriodogramWindowed<T>::Execute()
  100. {
  101.    int i,is;
  102. #ifdef _DEBUG
  103.    *DebugFile 
  104.       << "In BartlettPeriodogramWindowed::Execute" 
  105.       << endl;
  106. #endif
  107.    if(Processing_Completed) return(_MES_AOK);
  108.    if(PassNumber < Hold_Off) return (_MES_AOK);
  109.    //--------------------------------
  110.    //  Get pointers for buffers
  111.    T *in_sig_ptr = GET_INPUT_PTR(In_Sig);
  112.    int samps_avail = Block_Size;
  113.    while(Samps_Needed <= samps_avail)
  114.    {
  115.       //  The new input block has enough samples
  116.       //  to finish a segment.
  117.       //  Fill up FFT buffer by getting 
  118.       //  Samps_Needed input samples.
  119.       for(is=Samps_Needed; is>0; is--){
  120.          Time_Seg[Seg_Len - is] = *in_sig_ptr++;
  121.       }
  122.       samps_avail -= Samps_Needed;
  123.       if(Using_Window){
  124.          for(is=0; is<Seg_Len; is++){
  125.             Win_Time_Seg[is] = 
  126.                   Window_Taps[is]*Time_Seg[is];
  127.          }
  128.       }
  129.       else{
  130.          for(is=0; is<Seg_Len; is++){
  131.             Win_Time_Seg[is] = Time_Seg[is];
  132.          }
  133.       }
  134.       //  Perform FFT
  135.       FFT<double>(   Win_Time_Seg,
  136.                      Freq_Seg,
  137.                      Seg_Len,
  138.                      Fft_Len);
  139.       for(i=0; i<Seg_Len; i++){
  140.          Sample_Spectrum[i] += Samp_Intvl * 
  141.                   std::norm(Freq_Seg[i])/Seg_Len;
  142.       }
  143.       for(is=0; is<Seg_Len; is++){
  144.          Time_Seg[is] = 0.0;
  145.       }
  146.       Samps_Needed = Seg_Len;
  147.       Segs_In_Est++;
  148.       // is it time to dump the results?
  149.       if(Segs_In_Est == Num_Segs_To_Avg){
  150.          for(i=0; i<Seg_Len; i++){
  151.             Sample_Spectrum[i] /= 
  152.                         double(Num_Segs_To_Avg);
  153.          }
  154.          DumpSpectrum(  Sample_Spectrum,
  155.                         Fft_Len,
  156.                         Delta_F,
  157.                         Freq_Norm_Factor,
  158.                         Output_In_Decibels,
  159.                         Plot_Two_Sided,
  160.                         Psd_File);
  161.          Processing_Completed = true;
  162.          Psd_File->close();
  163.          if(Halt_When_Completed){
  164. #ifdef _DEBUG
  165.             *DebugFile << "Execution halted by " 
  166.                        << GetModelName() << endl;
  167. #endif
  168.             exit(0);
  169.          }
  170.       }
  171.    }// end of while
  172.    //---------------------------------------------------
  173.    //  The number of avail new samples is not sufficient
  174.    //  to finish a segment.  Copy the avaialble samples 
  175.    //  and then wait for the next pass to get some more.
  176.    for(is=0; is<samps_avail; is++){
  177.       Time_Seg[Seg_Len - Samps_Needed + is] = 
  178.                                        *in_sig_ptr++;
  179.    }
  180.    Samps_Needed -= samps_avail;
  181.    return(_MES_AOK);
  182. }
  183. //======================================================
  184. template BartlettPeriodogramWindowed<float>;