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

void paz::setup ()
{
  if (!killfile.empty ())
  {
    FILE *fptr = fopen (killfile.c_str (), "r");

    if (!fptr)
      throw Error (FailedSys, "paz::setup", "fopen " + killfile);

    char *buffer = new char[4096];
    while (fgets (buffer, 4096, fptr))
    {
      char* key = strtok (buffer, whitespace);
      while (key)
      {
	chans_to_zero.push_back (fromstring<unsigned>(key));
	key = strtok (NULL, whitespace);
      }
    }
  }

  zero_channels = chans_to_zero.size()
    || freqs_to_zero.size() || freq_ranges_to_zero.size();

  bool zero_subints = subs_to_zero.size ();

  if (zero_intersection && !(zero_channels && zero_subints))
    throw Error (InvalidState, "paz::setup",
		 "must use -I with both -s|S|k *and* -w|W");
  
  if (mow_all_subints || subints_to_mow.size())
    mower = new Pulsar::LawnMower;
  
  if (median_zap || median_zap_bybin || median_zap_window)
  {
    median_zapper = new Pulsar::ChannelZapMedian;
    median_zapper->set_bybin (median_zap_bybin);
    median_zapper->set_paz_report (true);

    if (median_zap_window)
      median_zapper->set_window_size (median_zap_window);
  }

  if (edge_zap_percent < 0.0 || edge_zap_percent >= 100.0)
    throw Error (InvalidState, "paz::setup",
		 "invalid parameter to option -E");

  if (!std_filename.empty())
  {
    Reference::To<Pulsar::Archive> data = Pulsar::Archive::load (std_filename);
    data->pscrunch ();
    data->fscrunch ();
    data->tscrunch ();
    thestd = data->get_Profile (0, 0, 0);
    
    standard_snr.set_standard (thestd);
    Pulsar::Profile::snr_strategy.get_value().set (&standard_snr,
						   &Pulsar::StandardSNR::get_snr);
  }
}

void psrpca::setup ()
{
  arrival->set_shift_estimator ( new PhaseGradShift );
  try
  {
    std_archive->fscrunch();
    std_archive->pscrunch();
    arrival->set_standard ( std_archive );
  } catch (Error& error)
  {
    cerr << error << endl;
    cerr << "psrpca::setup setting standard failed" << endl;
    exit ( -1 );
  }
}

/**
 * Zap a specified range of bins by substituting their values with the
 * baseline mean for each profile.
 */
void zap_bins(Pulsar::Archive* arch, Pulsar::Archive* old_arch,
    const unsigned s, const unsigned start_bin, const unsigned end_bin)
{
  const unsigned npol  = old_arch->get_npol();
  const unsigned nchan = old_arch->get_nchan();

  // Get the baseline region from what is displayed on screen (pscrunched
  // and fscrunched integration).
  Reference::To<Pulsar::Archive> copy = old_arch->clone();
  copy->pscrunch();
  copy->fscrunch();

  Reference::To<Pulsar::PhaseWeight> baseline_weights =
    copy->get_Profile(s,0,0)->baseline();

  vector<float> weights;
  baseline_weights->get_weights(weights);

  for (unsigned ipol = 0; ipol < npol; ++ipol) {
    for (unsigned ichan = 0; ichan < nchan; ++ichan) {
      // Apply the baseline region to each profile and get the mean.
      baseline_weights->set_Profile(old_arch->get_Profile(s,ipol,ichan));
      const double baseline_mean = baseline_weights->get_mean().get_value();
      float* bins = old_arch->get_Profile(s,ipol,ichan)->get_amps();

      // Set all the bins between start_bin and end_bin to the mean of the
      // baseline.
      for (unsigned ibin = start_bin; ibin <= end_bin; ++ibin) {
        bins[ibin] = baseline_mean;
      }
    }
  }

  *arch = *old_arch;
  arch->set_dispersion_measure(0);
  arch->pscrunch();
  arch->fscrunch();
  arch->remove_baseline();
}

//.........这里部分代码省略.........
	break;
      case 'E':

	try {
	  new_eph = factory<Pulsar::Parameters> (optarg);
	}
	catch (Error& error) {
	  cerr << "Could not load new ephemeris from " << optarg << endl;
	  return -1;
	}

	command += " -E";
	break;
      case 'T':
	tscr = true;
	command += " -T";
	break;
      case 'F':
	fscr = true;
	command += " -F";
	break;
      case 'p':
	pscr = true;
	command += " -p";
	break;
      case 'I':
	invint = true;
	pscr = false;
	command += " -I";
	break;
      case 'f':
	fscr = true;
	if (sscanf(optarg, "%d", &fscr_fac) != 1) {
	  cout << "That is not a valid fscrunch factor" << endl;
	  return -1;
	}
	command += " -f ";
	command += optarg;
	break;
	
      case 'n':

	reflections.add_reflection( optarg[0] );

	command += " -n ";
	command += optarg;
	break;

      case 'o':
	new_cfreq = true;
	if (sscanf(optarg, "%lf", &new_fr) != 1) {
	  cout << "That is not a valid centre frequency" << endl;
	  return -1;
	}
	command += " -o ";
	command += optarg;
	break;
      case 't':
	tscr = true;
	if (sscanf(optarg, "%d", &tscr_fac) != 1) {
	  cout << "That is not a valid tscrunch factor" << endl;
	  return -1;
	}
	command += " -t ";
	command += optarg;
	break;

//.........这里部分代码省略.........
      centre_model = true;
      break;

    case 'p':
      rotate_peak = true;
      break;

    case 'R':
      rotate_amount = atof(optarg);
      break;

    case 'a':
      align = true;
      break;

    case 'j':
      details_filename = optarg;
      break;

   default:
      cerr << "invalid param '" << c << "'" << endl;
    }

  if (!pgdev.empty())
  {
    cpgopen(pgdev.c_str());
    cpgask(0);
    cpgsvp(0.1, 0.9, 0.1, 0.9);
  }
  
  Reference::To<Archive> archive = Archive::load (argv[optind]);

  // preprocess
  archive->fscrunch();
  archive->tscrunch();
  archive->pscrunch();

  // phase up as requested
  if (centre_model)
    archive->centre();
  else if (rotate_peak)
    archive->centre_max_bin(0.0);
  
  if (rotate_amount != 0.0)
    archive->rotate_phase(-rotate_amount);
  
  archive->remove_baseline();

  // load from file if specified
  if (!model_filename_in.empty())
  {
    cerr << "paas: loading model from " << model_filename_in << endl;
    model.load(model_filename_in.c_str());
    
    // align to profile first if asked
    if (align)
      model.align(archive->get_Integration(0)->get_Profile(0,0));
  }

  // add any new components specified
  unsigned ic, nc=new_components.size();
  double centre, concentration, height;
  int name_start;
  for (ic=0; ic < nc; ic++)
  {
    if (sscanf(new_components[ic].c_str(), 

int main(int argc, char* argv[]) try
{
  if (argc < 2) {
    usage();
    return EXIT_SUCCESS;
  }
  
  int gotc = 0;
  while ((gotc = getopt(argc, argv, "hvV")) != -1) {
    switch (gotc) {
    case 'h':
      usage();
      return EXIT_SUCCESS;
    case 'V':
      Pulsar::Archive::set_verbosity(3);
      break;
    case 'v':
      Pulsar::Archive::set_verbosity(2);
      break;
    }
  }
  
  if (optind >= argc)
    {
      cerr << "pazi: please specify filename" << endl;
      return -1;
    }
  
  string filename = argv[optind];
  string extension = filename.substr(filename.length() - 2, 2);
  
  if (extension == "rf")
    extension = "rz";
  else if (extension == "cf")
    extension = "cz";
  else
    extension = "pazi";
  
  string write_filename = filename + ".";
  write_filename += extension;
  
  cerr << "pazi: loading data" << endl;
  base_archive = Archive::load(filename);
  
  if (base_archive->get_npol() == 4)
  {
    original_state = base_archive->get_state();
    base_archive->convert_state( Signal::Stokes );
  }
  
  backup_archive = base_archive->clone();
  
  cerr << "pazi: making fscrunched clone" << endl;
  mod_archive = base_archive->clone();
  mod_archive->pscrunch();
  mod_archive->remove_baseline();
  mod_archive->dedisperse();
  mod_archive->fscrunch();

  scrunched_archive = mod_archive->clone();
  scrunched_archive->tscrunch();

  ranges.second = get_max_value(base_archive, plot_type);
  positive_direction = base_archive->get_bandwidth() < 0.0;

  time_orig_plot = factory.construct("time");
  time_mod_plot = factory.construct("time");
  time_fui = time_mod_plot->get_frame_interface();

  freq_orig_plot = factory.construct("freq");
  freq_mod_plot = factory.construct("freq");
  freq_fui = freq_mod_plot->get_frame_interface();

  total_plot = factory.construct("flux");
  total_plot->configure("info=1");

  subint_orig_plot = new ProfilePlot;
  subint_mod_plot = new ProfilePlot;
  subint_fui = subint_mod_plot->get_frame_interface();

  subint_orig_plot->configure("info=1");
  subint_mod_plot->configure("info=1");

  unsigned window = 0;
  char device [8];

  for (unsigned i=0; i<2; i++) do {
    window ++;
    snprintf (device, 8, "%u/XS", window);
  }
  while ( cpgopen (device) < 0 );

  cpgask(0);

  cerr << endl << "Total S/N = " <<
    scrunched_archive->get_Profile(0,0,0)->snr() << endl << endl;

  total_plot->plot(scrunched_archive);
  cpgslct(1);
  time_orig_plot->plot(mod_archive);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    {
      for (unsigned ic=0; ic<nchan; ic++)
      {
	FreqRange chan;
	chan.lo=arch->get_Integration(0)->get_centre_frequency(ic) 
	  - chan_bw/2.0;
	chan.hi=arch->get_Integration(0)->get_centre_frequency(ic) 
	  + chan_bw/2.0;
	if (chan.lo > freq_ranges_to_zero[i].lo 
	    && chan.lo < freq_ranges_to_zero[i].hi)
	  mask[ic] = 0.0;
	else if (chan.hi > freq_ranges_to_zero[i].lo 
		 && chan.hi < freq_ranges_to_zero[i].hi)
	  mask[ic] = 0.0;
	else if (freq_ranges_to_zero[i].lo > chan.lo 
		 && freq_ranges_to_zero[i].lo < chan.hi)
	  mask[ic] = 0.0;
      }
    }
    
    if (zero_intersection)
      zapper->zap_very_specific (arch, mask, subs_to_zero);
    else
      zapper->zap_specific (arch, mask);
  }
  
  if (edge_zap_percent)
  {
    float fraction = edge_zap_percent / 100.0;
    unsigned buffer = unsigned (float (nchan) * fraction);

    vector < float >mask (nchan, 0.0);

    for (unsigned i = buffer; i < (nchan - buffer); i++)
      mask[i] = 1.0;

    zapper->zap_specific (arch, mask);
  }

  if (ston_cutoff > 0.0)
  {
    double theston = 0.0;
    Reference::To<Pulsar::Archive> cloned = arch->clone();
    cloned->pscrunch ();
    for (unsigned isub = 0; isub < arch->get_nsubint (); isub++) {
      for (unsigned ichan = 0; ichan < arch->get_nchan (); ichan++) {
	theston = cloned->get_Profile (isub, 0, ichan)->snr ();
	if (theston < ston_cutoff) {
	  arch->get_Integration (isub)->set_weight (ichan, 0.0);
	}
      }
    }
  }

  if (dropout_sigma)
  {
    Reference::To < Pulsar::Archive > cloned = arch->clone ();
    Reference::To < Pulsar::Profile > testprof = 0;

    cloned->pscrunch ();
    cloned->fscrunch ();
    cloned->remove_baseline ();

    vector < double >mins;
    for (unsigned isub = 0; isub < arch->get_nsubint (); isub++) {
      testprof = cloned->get_Profile (isub, 0, 0);
      mins.push_back (fabs (testprof->min ()));
    }

    cloned->tscrunch ();

    double mean, vari, varm;

    testprof = cloned->get_Profile (0, 0, 0);
    testprof->stats (&mean, &vari, &varm);

    for (unsigned isub = 0; isub < arch->get_nsubint (); isub++) {
      if (mins[isub] > dropout_sigma * sqrt (vari)) {
	cerr << "Zapping integration " << isub << endl;
	arch->get_Integration (isub)->uniform_weight (0.0);
      }
    }
  }

  if (mower)
  {
    if (mow_all_subints)
    {
      subints_to_mow.resize (arch->get_nsubint ());
      for (unsigned isub = 0; isub < arch->get_nsubint (); isub++)
	subints_to_mow[isub] = isub;
    }

    for (unsigned isub = 0; isub < subints_to_mow.size (); isub++)
    {
      cerr << "paz: mowing subint " << subints_to_mow[isub] << endl;
      mower->transform (arch->get_Integration (subints_to_mow[isub]));
    }
  }
}

int main (int argc, char *argv[])
{
  bool  verbose    = false;
  bool  gaussian   = false;
  float factor     = 5.0;
  bool  log        = false;
  bool  shownoise  = false;
  int   method     = 0;
  int   bins       = 30;
  float cphs       = 0.0;
  float dcyc       = 0.0;
  float norm       = 0.0;

  int gotc = 0;

  while ((gotc = getopt(argc, argv, "qvVhlg:p:w:n:b:Gs")) != -1) {
    switch (gotc) {
    case 'h':
      usage ();
      return 0;

    case 'v':
      Pulsar::Archive::set_verbosity(2);
      verbose = true;
      break;

    case 'V':
      verbose = true;
      Pulsar::Archive::set_verbosity(3);
      break;

    case 'l':
      log = true;
      break;

    case 'g':
      factor = atof(optarg);
      break;

    case 'n':
      norm = atof(optarg);
      break;
      
    case 'p':
      cphs = atof(optarg);
      if (cphs > 0.0 && cphs < 1.0) {
	method = 2;
      }
      else {
	method = 1;
      }
      break;
      
    case 'w':
      dcyc = atof(optarg);
      break;

    case 'b':
      bins = atoi(optarg);
      break;

    case 'G':
      gaussian = true;
      break;

    case 's':
      shownoise = true;
      break;
      
    default:
      cout << "Unrecognised option: " << gotc << endl;
    }
  }
 
  vector<string> archives;
  vector<pulse>  pulses;
  vector<pulse>  noise;
  vector<pulse>  pgiants;
  vector<pulse>  ngiants;
  vector<hbin>   pdata;
  vector<hbin>   ndata;

  Reference::To<Pulsar::Archive> arch;

  for (int ai = optind; ai < argc; ai++)             // extract filenames
    dirglob (&archives, argv[ai]);
  
  try {
    
    cout << "Loading archives and extracting data..." << endl;
    
    for (unsigned i = 0; i < archives.size(); i++) {
      
      if (verbose) cerr << "Loading " << archives[i] << endl;

      arch = Pulsar::Archive::load(archives[i]);
      
      arch->fscrunch();        // remove frequency and polarisation resolution
      arch->pscrunch();
      arch->remove_baseline(); // Remove the baseline level
//.........这里部分代码省略.........