refclock_wwv.c

网络时间协议NTP 源码 版本v4.2.0b 该源码用于linux平台下

 *
 * This routine implements a virtual station process used to acquire
 * minute sync and to mitigate among the ten frequency and station
 * combinations. During minute sync acquisition the process probes each
 * frequency in turn for the minute pulse from either station, which
 * involves searching through the entire minute of samples. After
 * finding a candidate, the process searches only the seconds before and
 * after the candidate for the signal and all other seconds for the
 * noise.
 *
 * Students of radar receiver technology will discover this algorithm
 * amounts to a range gate discriminator. The discriminator requires
 * that the peak minute pulse amplitude be at least 2000 and the SNR be
 * at least 6 dB. In addition after finding a candidate, The peak second
 * pulse amplitude must be at least 2000, the SNR at least 6 dB and the
 * difference between the current and previous epoch must be less than
 * 7.5 ms, which corresponds to a frequency error of 125 PPM. A compare
 * counter keeps track of the number of successive intervals which
 * satisfy these criteria.
 *
 * Note that, while the minute pulse is found by by the discriminator,
 * the actual value is determined from the second epoch. The assumption
 * is that the discriminator peak occurs about 800 ms into the second,
 * so the timing is retarted to the previous second epoch.
 */
static void
wwv_qrz(
	struct peer *peer,	/* peer structure pointer */
	struct sync *sp,	/* sync channel structure */
	int	pdelay		/* propagation delay (samples) */
	)
{
	struct refclockproc *pp;
	struct wwvunit *up;
	char tbuf[80];		/* monitor buffer */
	long epoch, fpoch;
	int isgood;

	pp = peer->procptr;
	up = (struct wwvunit *)pp->unitptr;

	/*
	 * Find the sample with peak energy, which defines the minute
	 * epoch. If a sample has been found with good energy,
	 * accumulate the noise energy for all except the second before
	 * and after that position.
	 *
	 * If the seconds sync lamp is lit, use that as the sample epoch
	 * within the second; otherwise, use the minutes sync peak.
	 * Little bit of class here.
	 */
	if (up->epomax > STHR && up->eposnr > SSNR) {
		fpoch = up->mphase % SECOND - up->tepoch;
		if (fpoch < 0)
			fpoch += SECOND;
	} else {
		fpoch = pdelay + SYNSIZ;
	}
	epoch = up->mphase - fpoch;
	if (epoch < 0)
		epoch += MINUTE;
	if (sp->amp > sp->maxeng) {
		sp->maxeng = sp->amp;
		sp->pos = epoch;
	}
	if (abs((epoch - sp->lastpos) % MINUTE) > SECOND)
		sp->noieng += sp->amp;

	/*
	 * At the end of the minute, determine the epoch of the
	 * sync pulse, as well as the SNR and difference between
	 * the current and previous epoch, which represents the
	 * intrinsic frequency error plus jitter.
	 */
	if (up->mphase == 0) {
		sp->synmax = sp->maxeng;
		sp->synsnr = wwv_snr(sp->synmax, sp->noieng / (MINUTE -
		    2. * SECOND));
		epoch = (sp->pos - sp->lastpos) % MINUTE;

		/*
		 * If not yet in minute sync, we have to do a little
		 * dance to find a valid minute sync pulse, emphasis
		 * valid.
		 */
		isgood = sp->synmax > ATHR && sp->synsnr > ASNR;
		switch (sp->count) {

		/*
		 * In state 0 the station was not heard during the
		 * previous probe. Look for the biggest blip greater
		 * than the amplitude threshold in the minute and assume
		 * that the minute sync pulse. We're fishing here, since
		 * the range gate has not yet been determined. If found,
		 * bump to state 1.
		 */
		case 0:
			if (sp->synmax >= ATHR)
				sp->count++;
			break;

		/*
		 * In state 1 a candidate blip has been found and the
		 * next minute has been searched for another blip. If
		 * none are found acceptable, drop back to state 0 and
		 * hunt some more. Otherwise, a legitimate minute pulse
		 * may have been found, so bump to state 2.
		 */
		case 1:
			if (!isgood) {
				sp->count = 0;
				break;
			}
			sp->count++;
			break;

		/*
		 * In states 2 and above, continue to groom samples as
		 * before and drop back to state 0 if the groom fails.
		 * If it succeeds, set the epoch and bump to the next
		 * state until reaching the threshold, if ever.
		 */
		default:
			if (!isgood || abs(epoch) > AWND * MS) {
				sp->count = 0;
				break;
			}
			sp->mepoch = sp->pos;
			sp->count++;
			break;
		}
		sp->metric = wwv_metric(sp);
		if (pp->sloppyclockflag & CLK_FLAG4) {
			sprintf(tbuf,
			    "wwv8 %d %3d %s %d %5.0f %5.1f %5.0f %5ld %5d %ld",
			    up->port, up->gain, sp->refid, sp->count,
			    sp->synmax, sp->synsnr, sp->metric, sp->pos,
			    up->tepoch, epoch);
			record_clock_stats(&peer->srcadr, tbuf);
#ifdef DEBUG
			if (debug)
				printf("%s\n", tbuf);
#endif
		}
		sp->lastpos = sp->pos;
		sp->maxeng = sp->noieng = 0;
	}
}


/*
 * wwv_endpoc - identify and acquire second sync pulse
 *
 * This routine is called at the end of the second sync interval. It
 * determines the second sync epoch position within the interval and
 * disciplines the sample clock using a frequency-lock loop (FLL).
 *
 * Second sync is determined in the RF input routine as the maximum
 * over all 8000 samples in the second comb filter. To assure accurate
 * and reliable time and frequency discipline, this routine performs a
 * great deal of heavy-handed heuristic data filtering and grooming.
 */
static void
wwv_endpoc(
	struct peer *peer,	/* peer structure pointer */
	int epopos		/* epoch max position */
	)
{
	struct refclockproc *pp;
	struct wwvunit *up;
	static int epoch_mf[3]; /* epoch median filter */
 	static int xepoch;	/* last second epoch */
 	static int zepoch;	/* last averaging interval epoch */
	static int syncnt;	/* run length counter */
	static int maxrun;	/* longest run length */
	static int mepoch;	/* longest run epoch */
	static int avgcnt;	/* averaging interval counter */
	static int avginc;	/* averaging ratchet */
	static int iniflg;	/* initialization flag */
	char tbuf[80];		/* monitor buffer */
	double dtemp;
	int tmp2;

	pp = peer->procptr;
	up = (struct wwvunit *)pp->unitptr;
	if (!iniflg) {
		iniflg = 1;
		memset((char *)epoch_mf, 0, sizeof(epoch_mf));
	}

	/*
	 * A three-stage median filter is used to help denoise the
	 * second sync pulse. The median sample becomes the candidate
	 * epoch.
	 */
	epoch_mf[2] = epoch_mf[1];
	epoch_mf[1] = epoch_mf[0];
	epoch_mf[0] = epopos;
	if (epoch_mf[0] > epoch_mf[1]) {
		if (epoch_mf[1] > epoch_mf[2])
			up->tepoch = epoch_mf[1];	/* 0 1 2 */
		else if (epoch_mf[2] > epoch_mf[0])
			up->tepoch = epoch_mf[0];	/* 2 0 1 */
		else
			up->tepoch = epoch_mf[2];	/* 0 2 1 */
	} else {
		if (epoch_mf[1] < epoch_mf[2])
			up->tepoch = epoch_mf[1];	/* 2 1 0 */
		else if (epoch_mf[2] < epoch_mf[0])
			up->tepoch = epoch_mf[0];	/* 1 0 2 */
		else
			up->tepoch = epoch_mf[2];	/* 1 2 0 */
	}

	/*
	 * If the signal amplitude or SNR fall below thresholds, dim the
	 * second sync lamp and start over. If no stations are heard we
	 * are either in a probe cycle or the ions are dim. In that case
	 * allow the amplitude and SNR to discharge and go no further. 
	 */
	if (up->epomax < STHR || up->eposnr < SSNR) {
		up->status &= ~(SSYNC | FGATE);
		avgcnt = syncnt = maxrun = 0;
		return;
	}
	if (!(up->status & (SELV | SELH)))
		return;

	avgcnt++;

	/*
	 * If the epoch candidate is the same as the last one, increment
	 * the compare counter. If not, save the length and epoch of the
	 * current run for use later and reset the counter.
	 */
	tmp2 = (up->tepoch - xepoch) % SECOND;
	if (tmp2 == 0) {
		syncnt++;
	} else {
		if (maxrun > 0 && mepoch == xepoch) {
			maxrun += syncnt;
		} else if (syncnt > maxrun) {
			maxrun = syncnt;
			mepoch = xepoch;
		}
		syncnt = 0;
	}
	if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status & (SSYNC |
	    MSYNC))) {
		sprintf(tbuf,
		    "wwv1 %04x %5.0f %5.1f %5d %5d %4d %4d",
		    up->status, up->epomax, up->eposnr, up->tepoch,
		    tmp2, avgcnt, syncnt);
		record_clock_stats(&peer->srcadr, tbuf);
#ifdef DEBUG
		if (debug)
			printf("%s\n", tbuf);
#endif /* DEBUG */
	}

	/*
	 * The sample clock frequency is disciplined using a first order
	 * feedback loop with time constant consistent with the Allan
	 * intercept of typical computer clocks.
	 *
	 * The frequency update is calculated from the epoch change in
	 * 125-us units divided by the averaging interval in seconds.
	 * The averaging interval affects other receiver functions,
	 * including the the 1000/1200-Hz comb filter and codec clock
	 * loop. It also affects the 100-Hz subcarrier loop and the bit
	 * and digit comparison counter thresholds.
	 */
	if (avgcnt < up->avgint) {
		xepoch = up->tepoch;
		return;
	}

	/*
	 * During the averaging interval the longest run of identical
	 * epoches is determined. If the longest run is at least 10
	 * seconds, the SSYNC bit is lit and the value becomes the
	 * reference epoch for the next interval. If not, the second
	 * sync lamp is dark and flashers set.
	 */
	if (maxrun > 0 && mepoch == xepoch) {
		maxrun += syncnt;
	} else if (syncnt > maxrun) {
		maxrun = syncnt;
		mepoch = xepoch;
	}
	xepoch = up->tepoch;
	if (maxrun > SCMP) {
		up->status |= SSYNC;
		up->yepoch = mepoch;
	} else {
		up->status &= ~SSYNC;
	}

	/*
	 * If the epoch change over the averaging interval is less than
	 * 1 ms, the frequency is adjusted, but clamped at +-125 PPM. If
	 * greater than 1 ms, the counter is decremented. If the epoch
	 * change is less than 0.5 ms, the counter is incremented. If
	 * the counter increments to +3, the averaging interval is
	 * doubled and the counter set to zero; if it decrements to -3,
	 * the interval is halved and the counter set to zero.
	 */
	if (maxrun == 0)
		mepoch = up->tepoch;
	dtemp = (mepoch - zepoch) % SECOND;
	if (up->status & FGATE) {
		if (abs(dtemp) < MAXFREQ * MINAVG) {
			up->freq += dtemp / (avgcnt * FCONST);
			if (up->freq > MAXFREQ)
				up->freq = MAXFREQ;
			else if (up->freq < -MAXFREQ)
				up->freq = -MAXFREQ;
			if (abs(dtemp) < MAXFREQ * MINAVG / 2.) {
				if (avginc < 3) {
					avginc++;
				} else {
					if (up->avgint < MAXAVG) {
						up->avgint <<= 1;
						avginc = 0;
					}
				}
			}
		} else {
			if (avginc > -3) {
				avginc--;
			} else {
				if (up->avgint > MINAVG) {
					up->avgint >>= 1;
					avginc = 0;
				}
			}
		}
	}
	if (pp->sloppyclockflag & CLK_FLAG4) {
		sprintf(tbuf,
		    "wwv2 %04x %4.0f %4d %4d %2d %4d %4.0f %7.2f",
		    up->status, up->epomax, mepoch, maxrun, avginc,
		    avgcnt, dtemp, up->freq * 1e6 / SECOND);
		record_clock_stats(&peer->srcadr, tbuf);
#ifdef DEBUG
		if (debug)
			printf("%s\n", tbuf);
#endif /* DEBUG */
	}
	up->status |= FGATE;
	zepoch = mepoch;
	avgcnt = syncnt = maxrun = 0;
}


/*
 * wwv_epoch - epoch scanner
 *
 * This routine extracts data signals from the 100-Hz subcarrier. It
 * scans the receiver second epoch to determine the signal amplitudes
 * and pulse timings. Receiver synchronization is determined by the
 * minute sync pulse detected in the wwv_rf() routine and the second
 * sync pulse detected in the wwv_epoch() routine. The transmitted
 * signals are delayed by the propagation delay, receiver delay and
 * filter delay of this program. Delay corrections are introduced
 * separately for WWV and WWVH. 
 *
 * Most communications radios use a highpass filter in the audio stages,
 * which can do nasty things to the subcarrier phase relative to the
 * sync pulses. Therefore, the data subcarrier reference phase is
 * disciplined using the hardlimited quadrature-phase signal sampled at
 * the same time as the in-phase signal. The phase tracking loop uses
 * phase adjustments of plus-minus one sample (125 us). Since this might
 * result in a phase slip of one cycle, the matched filter peak is the
 * maximum at the epoch and one cycle ahead or behind it. 
 */
static void
wwv_epoch(
	struct peer *peer	/* peer structure pointer */
	)
{
	struct refclockproc *pp;
	struct wwvunit *up;
	struct chan *cp;
	static double sigmin, sigzer, sigone, engmax, engmin;

	pp = peer->procptr;
	up = (struct wwvunit *)pp->unitptr;

	/*
	 * Sample the minute sync pulse energy at epoch 800 for both the
	 * WWV and WWVH stations. This will be used later for channel
	 * and station mitigation. Note that the seconds epoch is set
	 * here well before the end of the second to make sure we never
	 * set the epoch backwards.
	 */
	if (up->rphase == 800 * MS) {
		up->repoch = up->yepoch;
		cp = &up->mitig[up->achan];
		cp->wwv.syneng = cp->wwv.amp;
		cp->wwvh.syneng = cp->wwvh.amp;
	}

	/*
	 * Sample the I and Q data channels at epoch 200 ms. Use the
	 * signal energy as the peak to compute the SNR. Use the Q
	 * sample to adjust the 100-Hz reference oscillator phase.
	 */
	if (up->rphase == 200 * MS) {
		engmax = sqrt(up->irig * up->irig + up->qrig *
		    up->qrig);
		up->datpha = up->qrig / up->avgint;
		if (up->datpha >= 0) {
			up->datapt++;
			if (up->datapt >= 80)
				up->datapt -= 80;
		} else {
			up->datapt--;
			if (up->datapt < 0)
				up->datapt += 80;
		}
	}

	/*
	 * Use the signal amplitude at epoch 15 ms as the noise floor.
	 * This give a guard time of +-15 ms from the beginning of the
	 * second until the pulse rises at 30 ms. There is a compromise
	 * here; we want to delay the sample as long as possible to give
	 * the radio time to change frequency and the AGC to stabilize,