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

/* 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);
}

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();
}

//.........这里部分代码省略.........
      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

      scan_pulses(arch, pulses, method, cphs, dcyc);
      if (shownoise)
	scan_pulses(arch, noise, 2, arch->find_min_phase(dcyc), dcyc);
    }
    
    float offset = 2.0;
    
    if (shownoise) {
      offset = noise[0].flx;
      for (unsigned i = 1; i < noise.size(); i++)
	if (noise[i].flx < offset) offset = noise[i].flx;
      offset = fabs(offset);
      offset /= norm;
    }
    
    float threshold = find_giants(pulses, pgiants, factor, norm, offset);
    cout << "Detection threshold is roughly " << threshold << endl;

    display_giants(pgiants);

    prob_hist(pulses, pdata, bins);

    float submin = pdata[0].x;
    float submax = pdata[0].x;
    
    for (unsigned i = 1; i < pdata.size(); i++) {
      if (pdata[i].x > submax) submax = pdata[i].x;
      if (pdata[i].x < submin) submin = pdata[i].x;
    }

    if (shownoise) {
      find_giants(noise, ngiants, factor, norm, offset);
      prob_hist(noise, ndata, bins, submin, submax);
    }