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

3G开发

开发平台:

Visual C++

  1. //
  2. //  File = yule_walk_psd.cpp
  3. //
  4. #include <stdlib.h>
  5. #include "parmfile.h"
  6. #include "model_graph.h"
  7. #include "yule_walk_psd.h"
  8. #include "trianglr.h"
  9. #include "hamming.h"
  10. #include "hann.h"
  11. #include "fft_T.h"
  12. #include "dump_spect.h"
  13. #include "autometh.h"
  14. #include "levin.h"
  15. #include "ar_spec.h"
  16. #ifdef _DEBUG
  17.   extern ofstream *DebugFile;
  18. #endif
  19. #define _NO_ZERO_ENDS 0
  20. #define _ZERO_ENDS 1
  21. extern ParmFile *ParmInput;
  22. extern int PassNumber;
  23. //======================================================
  24. YuleWalkerPsdEstim::YuleWalkerPsdEstim(  
  25.                               char* instance_name,
  26.                               PracSimModel* outer_model,
  27.                               Signal<float>* in_sig )
  28.          :PracSimModel( instance_name,
  29.                         outer_model )
  30. {
  31.    MODEL_NAME(YuleWalkerPsdEstim);
  32.    OPEN_PARM_BLOCK;
  33.    GET_INT_PARM(Seg_Len);
  34.    GET_INT_PARM(Ar_Order);
  35.    GET_INT_PARM(Hold_Off);
  36.    GET_INT_PARM(Num_Freq_Pts);
  37.    Psd_File_Name = new char[64];
  38.    strcpy(Psd_File_Name, "");
  39.    GET_STRING_PARM(Psd_File_Name);
  40.    GET_DOUBLE_PARM(Freq_Norm_Factor);
  41.    GET_BOOL_PARM(Output_In_Decibels);
  42.    GET_BOOL_PARM(Plot_Two_Sided);
  43.    GET_BOOL_PARM(Halt_When_Completed);
  44.    In_Sig = in_sig;
  45.    MAKE_INPUT(In_Sig);
  46.    Time_Seg = new double[Seg_Len];
  47.    Sample_Spectrum = new double[Seg_Len];
  48.    for(int is=0; is<Seg_Len; is++){
  49.       Sample_Spectrum[is] = 0.0;
  50.    }
  51.    Processing_Completed = false;
  52. }
  53. //======================================================
  54. YuleWalkerPsdEstim::~YuleWalkerPsdEstim( void ){ };
  55. //======================================================
  56. void YuleWalkerPsdEstim::Initialize(void)
  57. {
  58.    Segs_In_Est = 0;
  59.    Samps_Needed = Seg_Len;
  60.    Block_Size = In_Sig->GetBlockSize();
  61.    Samp_Intvl = In_Sig->GetSampIntvl();
  62. };
  63. //======================================================
  64. int YuleWalkerPsdEstim::Execute()
  65. {
  66.    int is;
  67.    int error_status;
  68. #ifdef _DEBUG
  69.    *DebugFile << "In YuleWalkerPsdEstim::Execute" 
  70.               << endl;
  71. #endif
  72.    if(Processing_Completed) return(_MES_AOK);
  73.    if(PassNumber < Hold_Off) return (_MES_AOK);
  74.    //--------------------------------
  75.    //  Get pointers for buffers
  76.    float *in_sig_ptr = GET_INPUT_PTR(In_Sig);
  77.    int samps_avail = Block_Size;
  78.    while(Samps_Needed <= samps_avail){
  79.       //  The new input block has enough samples to
  80.       //  finish a segment.
  81.       //  Fill up FFT buffer by getting Samps_Needed
  82.       //  input samples.
  83.       for(is=Samps_Needed; is>0; is--){
  84.          Time_Seg[Seg_Len - is] = *in_sig_ptr++;
  85.       }
  86.       samps_avail -= Samps_Needed;
  87.       AutocorrMethCorrMtx *corr_mtx = 
  88.             new AutocorrMethCorrMtx(  
  89.                   Time_Seg, 
  90.                   Seg_Len, 
  91.                   Ar_Order);
  92.       double *correl_vec = corr_mtx->GetCol(1);
  93.       double *a_vec = new double[Ar_Order+1];
  94.       double driving_noise_var;
  95.       error_status = LevinsonRecursion(   
  96.          correl_vec,
  97.          Ar_Order,
  98.          a_vec,
  99.          &driving_noise_var);
  100.       ArSpectrum *ar_spectrum = new ArSpectrum( 
  101.                Ar_Order,
  102.                a_vec,
  103.                Samp_Intvl,
  104.                driving_noise_var, //true_ar_drv_var );
  105.                Num_Freq_Pts,
  106.                1.0/(2*Samp_Intvl));
  107.       ar_spectrum->DumpSpectrum( Psd_File_Name, true);
  108.       if(Halt_When_Completed)
  109.       {
  110. #ifdef _DEBUG
  111.          *DebugFile << "Execution halted by " 
  112.                     << GetModelName() << endl;
  113. #endif
  114.          exit(0);
  115.       }
  116.    }// end of while
  117.    //---------------------------------------------------
  118.    //  The number of avail new samples is not sufficient
  119.    //  to finish a segment.  Copy the avaialble samples
  120.    //  and then wait for the next pass to get some more.
  121.    for(is=0; is<samps_avail; is++){
  122.       Time_Seg[Seg_Len - Samps_Needed + is] = 
  123.                                        *in_sig_ptr++;
  124.    }
  125.    Samps_Needed -= samps_avail;
  126.    return(_MES_AOK);
  127. }
  128. //======================================================