C++ To::get_amps方法代码示例

/* This algorithm uses zero padding in the fourier domain to
   interpolate the cross correlation function in the time domain. The
   technique is described at:
  
   http://www.dspguru.com/howto/tech/zeropad2.htm 
*/
Estimate<double> Pulsar::ZeroPadShift::get_shift () const
{
  // First compute the standard cross correlation function:

  Reference::To<Pulsar::Profile> ptr = observation->clone();
  Reference::To<Pulsar::Profile> stp = standard->clone();

  // Remove the baseline
  *ptr -= ptr->mean(ptr->find_min_phase(0.15));

  // Perform the correlation
  ptr->correlate (standard);

  // Remove the baseline
  *ptr -= ptr->mean(ptr->find_min_phase(0.15));
  
  vector< Estimate<float> > correlation;

  const unsigned nbin = observation->get_nbin();
 
  for (unsigned i = 0; i < nbin; i++) {
    correlation.push_back(ptr->get_amps()[i]);
  }
  
  vector< Estimate<float> > interpolated;
  
  interpolated.resize(correlation.size() * interpolate);
  
  // Perform the zero-pad interpolation
  
  fft::interpolate(interpolated, correlation);
  
  // Find the peak of the correlation function
  
  float maxval = 0.0;
  float maxloc = 0.0;
  
  for (unsigned i = 0; i < interpolated.size(); i++) {
    if (interpolated[i].val > maxval) {
      maxval = interpolated[i].val;
      maxloc = float(i) / interpolate;
    }
  }
  
  // Error estimate (???)
  float ephase = 1.0 / float(nbin);
  
  double shift = double(maxloc) / double(nbin);
  
  if (shift < -0.5)
    shift += 1.0;
  else if (shift > 0.5)
    shift -= 1.0;
  
  return Estimate<double>(shift,ephase*ephase);
}

Pulsar::Profile* differentiate (const Pulsar::Profile* profile, unsigned off=1)
{
  Reference::To<Pulsar::Profile> difference = profile->clone();

  unsigned nbin = profile->get_nbin();
  const float* amps = profile->get_amps();
  float* damps = difference->get_amps();

  for (unsigned ibin=0; ibin < nbin; ibin++)
    damps[ibin] = amps[(ibin+off)%nbin] - amps[ibin];

  return difference.release();
}

void Pulsar::SquareWave::get_transitions (const Profile* profile,
					  vector<unsigned>& up,
					  vector<unsigned>& down)
{
  unsigned nbin = profile->get_nbin();
  unsigned scrunch = 1;

  Reference::To<Profile> clone;
  if (use_nbin && nbin > use_nbin) {
    clone = profile->clone();
    scrunch = nbin/use_nbin;
    clone->bscrunch(scrunch);
    profile = clone;
    nbin = profile->get_nbin();
  }

  unsigned offset = (unsigned) (risetime * nbin);

cerr << "nbin=" << nbin << " offset=" << offset << endl;

  // differentiate the profile
  Reference::To<Profile> difference = differentiate (profile, offset);

  float* amps = difference->get_amps();

  // find the phase window in which the mean is closest to zero
  BaselineWindow window;
  window.set_find_mean (0.0);
  window.get_smooth()->set_turns (0.2);
  float zero = window.find_phase (nbin, amps);

  // get the noise statistics of the zero mean region
  double mean = 0;
  double variance = 0;
  difference->stats (zero, &mean, &variance);

cerr << "mean=" << mean << " rms=" << sqrt(variance) << endl;

  // check that the mean is actually zero
  double rms = sqrt(variance);
  if (mean > rms)
    throw Error (InvalidState, "Pulsar::SquareWave::get_transitions",
		 "mean=%lf > rms=%lf", mean, rms);

  float cutoff = threshold * rms;

  find_transitions (nbin, amps, up, cutoff);
  find_transitions (nbin, amps, down, -cutoff);
}

void Pulsar::SpectrumPlot::get_spectra (const Archive* data)
{
  unsigned nchan = data->get_nchan();

  spectra.resize(1);
  spectra[0].resize(nchan);

  Reference::To<ProfileStats> stats;
  Reference::To<TextInterface::Parser> parser;

  if (!expression.empty())
  {
    stats = new ProfileStats;
    parser = stats->get_interface ();
  }

  Reference::To<const Integration> subint;
  subint = get_Integration (data, isubint);

  for (unsigned ichan=0; ichan<nchan; ichan++)
  {
    Reference::To<const Profile> profile;
    profile = get_Profile (subint, ipol, ichan);
    if (profile -> get_weight() == 0.0)
      spectra[0][ichan] = 0.0;
    else if (stats)
    {
      stats->set_Profile (profile);
      string value = process( parser, expression );
      spectra[0][ichan] = fromstring<float>( value );
    }
    else if (ibin.get_integrate())
      spectra[0][ichan] = profile->sum() / (float)profile->get_nbin();
    else
      spectra[0][ichan] = profile->get_amps()[ibin.get_value()];
  }
}

void psrpca::finalize ()
{
  arrival->set_observation ( total );
  arrival->get_toas(toas);

  if ( remove_std_baseline )
    std_archive -> remove_baseline ();
  Reference::To<Profile> std_prof = std_archive->get_Profile(0, 0, 0);
  float *s_amps = std_prof->get_amps ();
  const float nbin = std_prof->get_nbin ();

  double scale, offset, snr;

  if ( total_count < nbin )
    cerr << "WARNING: psrpca::finalize - not enough observations provided, "
	    "covariance matrix will not have full rank" << endl;

  //total->remove_baseline();

  gsl_matrix *profiles = gsl_matrix_alloc ( (unsigned)nbin, total_count );

  for (unsigned i_subint = 0; i_subint < total->get_nsubint(); i_subint++ )
  {
    Reference::To<Profile> prof = total->get_Profile ( i_subint, 0, 0 );
    if ( apply_shift )
      prof->rotate_phase ( toas[i_subint].get_phase_shift() );
    snr = prof->snr ();

    //calculate the scale
    float *p_amps = prof->get_amps ();
    scale = 0.0;

    for ( unsigned i_bin = 0; i_bin < nbin; i_bin++ )
    {
      scale += s_amps[i_bin] * p_amps[i_bin];
    }
    scale = (prof->get_nbin()* scale - prof->sum() * std_prof->sum()) /
	    (prof->get_nbin()* std_prof->sumsq() - std_prof->sum() * std_prof->sum());

    // calculate the baseline offset
    offset = (scale * std_prof->sum() - prof->sum()) / nbin;

    if ( prof_to_std )
    {
      //match the profile to standard and subtract the standard
      if ( apply_offset )
	prof->offset  ( offset );
      if ( apply_scale )
	prof->scale ( 1.0/scale );
      prof->diff ( std_prof );

      double* damps;
      damps = new double [ (unsigned)nbin ];
      transform( prof->get_amps(), prof->get_amps() + (unsigned)nbin, damps, CastToDouble() );

      gsl_vector_const_view view = gsl_vector_const_view_array( damps, nbin );
      gsl_matrix_set_col ( profiles, i_subint, &view.vector );
      t_cov->add_Profile ( prof, snr );
    }
    else
    {// prof_to_std is false
      Reference::To<Profile> diff = prof->clone ();
      diff->set_amps ( std_prof->get_amps () );
      if ( apply_offset ) 
	diff->offset( -offset );
      if ( apply_scale ) 
	diff->scale (scale);
      diff->diff ( prof );
      diff->scale (-1);

      double* damps;
      damps = new double [ (unsigned)nbin ];
      transform( diff->get_amps(), diff->get_amps() + (unsigned)nbin, damps, CastToDouble() );

      gsl_vector_const_view view = gsl_vector_const_view_array( damps, nbin );
      gsl_matrix_set_col ( profiles, i_subint, &view.vector );

      t_cov->add_Profile ( diff, snr );
      prof->set_amps ( diff->get_amps() );
    }
  }

  covariance = gsl_matrix_alloc ( (int) nbin, (int) nbin );
  t_cov->get_covariance_matrix_gsl ( covariance );

  // write the covariance matrix and difference profiles
  FILE *out;
  if ( save_covariance_matrix )
  {
    out = fopen ( (prefix+"_covariance.dat").c_str(), "w" );
    gsl_matrix_fprintf(out, covariance, "%g");
    fclose ( out );
  } // save covariance matrix

  if ( save_diffs ) 
    total->unload ( prefix+"_diffs.ar" );

  //solve the eigenproblem
  gsl_matrix_view m = gsl_matrix_submatrix ( covariance, 0, 0, (int)nbin, (int)nbin );
  gsl_vector *eval = gsl_vector_alloc ( (int)nbin );
//.........这里部分代码省略.........

int scan_pulses(Reference::To<Pulsar::Archive> arch, vector<pulse>& data, 
		int method, float cphs, float dcyc)
{
  /* Find and store the flux and phase of each pulse in the file
     to the data vector. */
  
  pulse newentry;

  Reference::To<Pulsar::Profile> prof;
  
  newentry.file = arch->get_filename();

  double nm, nv, vm;
  
  int nbin = arch->get_nbin();
  int bwid = int(float(nbin) * dcyc);
  int cbin = 0;

  for (unsigned i = 0; i < arch->get_nsubint(); i++) {

    newentry.intg = i;
    
    prof = arch->get_Profile(i, 0, 0);
    
    prof->stats(prof->find_min_phase(), &nm, &nv, &vm);
    newentry.err = sqrt(nv);

    switch (method) {
      
    case 0: // Method of total flux
      
      newentry.flx = prof->sum();
      newentry.phs = prof->find_max_phase();
      break;
      
    case 1: // Method of maximum amplitude

      if (dcyc == 0.0) {
	newentry.flx = prof->max();
	newentry.phs  = prof->find_max_phase();
      }
      else {
	cbin = int(prof->find_max_phase(dcyc) * float(nbin));
	newentry.flx = prof->sum(cbin - bwid/2, cbin + bwid/2);
	newentry.phs  = prof->find_max_phase(dcyc);
      }
      break;

    case 2: // User-defined phase centre
      
      cbin = int(float(nbin) * cphs);

      if (dcyc == 0.0) {
	newentry.flx = (prof->get_amps())[cbin];
	newentry.phs = cphs;
      }
      else {
	newentry.flx = prof->sum(cbin - bwid/2, cbin + bwid/2);
	newentry.phs = cphs;
      }
      
      break;
      
    default:
      cerr << "No phase selection method chosen!" << endl;
    }
    
    data.push_back(newentry);
  }
  
  return data.size();
}

void psrspa::create_histograms ( Reference::To<Archive> archive )
{
  if ( verbose )
    cerr << "psrspa::create_histograms entered" << endl;
  // ensure Stokes parameters if creating polarisation histograms
  if ( create_polar_degree || create_polar_angle )
  {
    archive->convert_state ( Signal::Stokes );
    if ( verbose )
      cerr << "psrspa::create_histograms converted state of the archive to Stokes" << endl;
  }

  // auxillary vectors
  //vector< Estimate<float > > aux_vec_f;
  vector< Estimate<double > > aux_vec_d;

  // Full Stokes profile
  Reference::To<PolnProfile> profile;
  // Polarized flux
  Reference::To<Profile> P;
  P = new Profile;
  // Total flux
  float *T_amps = new float [ nbin ];
  float *P_amps = new float [ nbin ];

  unsigned bin_min, bin_max;

  if ( verbose && max_bscrunch > 1 )
    cerr << "psrspa::create_histograms entering the bscrunch loop for the " << archive->get_filename () << " archive" << endl;
  // TODO Uh, this should be fixed up. The first idea was to enable the phase resolved histograms to be bscrunch-aware as well, but I think it doesn't make much sense so I decided to keep only the max flux bscrunch aware. This can be done much more neatly probably in such a case
  for ( current_bscrunch = 1 ; current_bscrunch <= max_bscrunch ; current_bscrunch *= 2 )
  {
    // in each passage, we bscrunch by a factor of 2
    if ( current_bscrunch > 1 )
    {
      if ( verbose )
	cerr << "psrspa::create_histograms bscrunching the archive " << archive->get_filename () << " by a factor of 2" << endl;
      archive->bscrunch ( 2 );
    }
    if ( verbose )
      cerr << "psrspa::create_histograms entering the loop through subints of the " << archive->get_filename () << " archive" << endl;
    // loop through subints
    for ( unsigned isub = 0; isub < archive->get_nsubint (); isub++ )
    {
      if ( verbose )
	cerr << "psrspa::create_histograms creating necessary profiles for subint " << isub << " of " << archive->get_filename () << endl;
      if ( create_polar_angle || create_polar_degree && current_bscrunch == 1 )
      {
	profile = archive->get_Integration(isub)->new_PolnProfile(0);
	if ( verbose )
	  cerr << "psrspa::create_histograms retrieved PolnProfile for subint " << isub << " of " << archive->get_filename () << endl;
      }
      if ( create_polar_angle && current_bscrunch == 1 )
      {
	profile->get_orientation ( aux_vec_d, 0 );
	if ( verbose )
	  cerr << "psrspa::create_histograms retrieved polarisation angle for subint " << isub << " of " << archive->get_filename () << endl;
      }
      if ( create_polar_degree || create_flux || find_max_amp_in_range )
      {
	stats.set_profile ( archive->get_Profile ( isub, 0, 0 ) );
	b_sigma = sqrt ( stats.get_baseline_variance ().get_value () );
	T_amps = archive->get_Profile ( isub, 0, 0 )->get_amps (); 
	if ( verbose )
	  cerr << "psrspa::create_histograms retrieved total flux amps for subint " << isub << " of " << archive->get_filename () << endl;
	if ( create_polar_degree && current_bscrunch == 1 )
	{
	  profile->get_polarized ( P );
	  if ( verbose )
	    cerr << "psrspa::create_histograms retrieved polarized flux profile for subint "  << isub << " of " << archive->get_filename () << endl;
	  P_amps = P->get_amps ();
	  if ( verbose )
	    cerr << "psrspa::create_histograms retrieved polarized flux amps for subint " << isub << " of " << archive->get_filename () << endl;
	}
      }

      if ( verbose )
	cerr << "psrspa::create_histograms looping through the provided phase ranges for subint " << isub << " of " << archive->get_filename () << endl;
      unsigned curr_hist = 0;
      // loop through phase ranges
      for ( unsigned irange = 0; irange < phase_range.size () ; irange++ )
      {
	if ( irange%2 == 0)
	{
	  bin_min = unsigned( floor ( phase_range[irange] * float(nbin / current_bscrunch ) + 0.5 ) );
	  if ( verbose ) 
	    cerr << "psrspa::create_histograms set minimal bin to " << bin_min << endl;
	}
	else
	{
	  bin_max = unsigned( floor ( phase_range[irange] * float(nbin / current_bscrunch ) + 0.5 ) );
	  if ( bin_max == nbin )
	    bin_max = nbin - 1 ;
	  if ( verbose ) 
	    cerr << "psrspa::create_histograms set maximum bin to " << bin_max << endl;
	  // loop through bins in the given phase range
	  for ( unsigned ibin = bin_min; ibin <= bin_max; ibin++ )
	  {
	    if ( create_polar_angle && current_bscrunch == 1 )
	    {
//.........这里部分代码省略.........