refclock_wwv.c

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

 * wwv_poll - called by the transmit procedure
 *
 * This routine keeps track of status. If no offset samples have been
 * processed during a poll interval, a timeout event is declared. If
 * errors have have occurred during the interval, they are reported as
 * well. Once the clock is set, it always appears reachable, unless
 * reset by watchcat timeout.
 */
static void
wwv_poll(
	int	unit,		/* instance number (not used) */
	struct peer *peer	/* peer structure pointer */
	)
{
	struct refclockproc *pp;
	struct wwvunit *up;

	pp = peer->procptr;
	up = (struct wwvunit *)pp->unitptr;
	if (pp->coderecv == pp->codeproc)
		up->errflg = CEVNT_TIMEOUT;
	if (up->errflg)
		refclock_report(peer, up->errflg);
	up->errflg = 0;
	pp->polls++;
}


/*
 * wwv_rf - process signals and demodulate to baseband
 *
 * This routine grooms and filters decompanded raw audio samples. The
 * output signals include the 100-Hz baseband data signal in quadrature
 * form, plus the epoch index of the second sync signal and the second
 * index of the minute sync signal.
 *
 * There are two 1-s ramps used by this program. Both count the 8000
 * logical clock samples spanning exactly one second. The epoch ramp
 * counts the samples starting at an arbitrary time. The rphase ramp
 * counts the samples starting at the 5-ms second sync pulse found
 * during the epoch ramp.
 *
 * There are two 1-m ramps used by this program. The mphase ramp counts
 * the 480,000 logical clock samples spanning exactly one minute and
 * starting at an arbitrary time. The rsec ramp counts the 60 seconds of
 * the minute starting at the 800-ms minute sync pulse found during the
 * mphase ramp. The rsec ramp drives the seconds state machine to
 * determine the bits and digits of the timecode. 
 *
 * Demodulation operations are based on three synthesized quadrature
 * sinusoids: 100 Hz for the data signal, 1000 Hz for the WWV sync
 * signal and 1200 Hz for the WWVH sync signal. These drive synchronous
 * matched filters for the data signal (170 ms at 100 Hz), WWV minute
 * sync signal (800 ms at 1000 Hz) and WWVH minute sync signal (800 ms
 * at 1200 Hz). Two additional matched filters are switched in
 * as required for the WWV second sync signal (5 ms at 1000 Hz) and
 * WWVH second sync signal (5 ms at 1200 Hz).
 */
static void
wwv_rf(
	struct peer *peer,	/* peerstructure pointer */
	double isig		/* input signal */
	)
{
	struct refclockproc *pp;
	struct wwvunit *up;
	struct sync *sp, *rp;

	static double lpf[5];	/* 150-Hz lpf delay line */
	double data;		/* lpf output */
	static double bpf[9];	/* 1000/1200-Hz bpf delay line */
	double syncx;		/* bpf output */
	static double mf[41];	/* 1000/1200-Hz mf delay line */
	double mfsync;		/* mf output */

	static int iptr;	/* data channel pointer */
	static double ibuf[DATSIZ]; /* data I channel delay line */
	static double qbuf[DATSIZ]; /* data Q channel delay line */

	static int jptr;	/* sync channel pointer */
	static int kptr;	/* tick channel pointer */

	static int csinptr;	/* wwv channel phase */
	static double cibuf[SYNSIZ]; /* wwv I channel delay line */
	static double cqbuf[SYNSIZ]; /* wwv Q channel delay line */
	static double ciamp;	/* wwv I channel amplitude */
	static double cqamp;	/* wwv Q channel amplitude */

	static double csibuf[TCKSIZ]; /* wwv I tick delay line */
	static double csqbuf[TCKSIZ]; /* wwv Q tick delay line */
	static double csiamp;	/* wwv I tick amplitude */
	static double csqamp;	/* wwv Q tick amplitude */

	static int hsinptr;	/* wwvh channels phase */
	static double hibuf[SYNSIZ]; /* wwvh I channel delay line */
	static double hqbuf[SYNSIZ]; /* wwvh Q channel delay line */
	static double hiamp;	/* wwvh I channel amplitude */
	static double hqamp;	/* wwvh Q channel amplitude */

	static double hsibuf[TCKSIZ]; /* wwvh I tick delay line */
	static double hsqbuf[TCKSIZ]; /* wwvh Q tick delay line */
	static double hsiamp;	/* wwvh I tick amplitude */
	static double hsqamp;	/* wwvh Q tick amplitude */

	static double epobuf[SECOND]; /* epoch sync comb filter */
	static double epomax;	/* epoch sync amplitude buffer */
	static int epopos;	/* epoch sync position buffer */

	static int iniflg;	/* initialization flag */
	int	epoch;		/* comb filter index */
	int	pdelay;		/* propagation delay (samples) */
	double	dtemp;
	int	i;

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

	if (!iniflg) {
		iniflg = 1;
		memset((char *)lpf, 0, sizeof(lpf));
		memset((char *)bpf, 0, sizeof(bpf));
		memset((char *)mf, 0, sizeof(mf));
		memset((char *)ibuf, 0, sizeof(ibuf));
		memset((char *)qbuf, 0, sizeof(qbuf));
		memset((char *)cibuf, 0, sizeof(cibuf));
		memset((char *)cqbuf, 0, sizeof(cqbuf));
		memset((char *)csibuf, 0, sizeof(csibuf));
		memset((char *)csqbuf, 0, sizeof(csqbuf));
		memset((char *)hibuf, 0, sizeof(hibuf));
		memset((char *)hqbuf, 0, sizeof(hqbuf));
		memset((char *)hsibuf, 0, sizeof(hsibuf));
		memset((char *)hsqbuf, 0, sizeof(hsqbuf));
		memset((char *)epobuf, 0, sizeof(epobuf));
	}

	/*
	 * Baseband data demodulation. The 100-Hz subcarrier is
	 * extracted using a 150-Hz IIR lowpass filter. This attenuates
	 * the 1000/1200-Hz sync signals, as well as the 440-Hz and
	 * 600-Hz tones and most of the noise and voice modulation
	 * components.
	 *
	 * The subcarrier is transmitted 10 dB down from the carrier.
	 * The DGAIN parameter can be adjusted for this and to
	 * compensate for the radio audio response at 100 Hz.
	 *
	 * Matlab IIR 4th-order IIR elliptic, 150 Hz lowpass, 0.2 dB
	 * passband ripple, -50 dB stopband ripple.
	 */
	data = (lpf[4] = lpf[3]) * 8.360961e-01;
	data += (lpf[3] = lpf[2]) * -3.481740e+00;
	data += (lpf[2] = lpf[1]) * 5.452988e+00;
	data += (lpf[1] = lpf[0]) * -3.807229e+00;
	lpf[0] = isig * DGAIN - data;
	data = lpf[0] * 3.281435e-03
	    + lpf[1] * -1.149947e-02
	    + lpf[2] * 1.654858e-02
	    + lpf[3] * -1.149947e-02
	    + lpf[4] * 3.281435e-03;

	/*
	 * The 100-Hz data signal is demodulated using a pair of
	 * quadrature multipliers, matched filters and a phase lock
	 * loop. The I and Q quadrature data signals are produced by
	 * multiplying the filtered signal by 100-Hz sine and cosine
	 * signals, respectively. The signals are processed by 170-ms
	 * synchronous matched filters to produce the amplitude and
	 * phase signals used by the demodulator. The signals are scaled
	 * to produce unit energy at the maximum value.
	 */
	i = up->datapt;
	up->datapt = (up->datapt + IN100) % 80;
	dtemp = sintab[i] * data / (MS / 2. * DATCYC);
	up->irig -= ibuf[iptr];
	ibuf[iptr] = dtemp;
	up->irig += dtemp;

	i = (i + 20) % 80;
	dtemp = sintab[i] * data / (MS / 2. * DATCYC);
	up->qrig -= qbuf[iptr];
	qbuf[iptr] = dtemp;
	up->qrig += dtemp;
	iptr = (iptr + 1) % DATSIZ;

	/*
	 * Baseband sync demodulation. The 1000/1200 sync signals are
	 * extracted using a 600-Hz IIR bandpass filter. This removes
	 * the 100-Hz data subcarrier, as well as the 440-Hz and 600-Hz
	 * tones and most of the noise and voice modulation components.
	 *
	 * Matlab 4th-order IIR elliptic, 800-1400 Hz bandpass, 0.2 dB
	 * passband ripple, -50 dB stopband ripple.
	 */
	syncx = (bpf[8] = bpf[7]) * 4.897278e-01;
	syncx += (bpf[7] = bpf[6]) * -2.765914e+00;
	syncx += (bpf[6] = bpf[5]) * 8.110921e+00;
	syncx += (bpf[5] = bpf[4]) * -1.517732e+01;
	syncx += (bpf[4] = bpf[3]) * 1.975197e+01;
	syncx += (bpf[3] = bpf[2]) * -1.814365e+01;
	syncx += (bpf[2] = bpf[1]) * 1.159783e+01;
	syncx += (bpf[1] = bpf[0]) * -4.735040e+00;
	bpf[0] = isig - syncx;
	syncx = bpf[0] * 8.203628e-03
	    + bpf[1] * -2.375732e-02
	    + bpf[2] * 3.353214e-02
	    + bpf[3] * -4.080258e-02
	    + bpf[4] * 4.605479e-02
	    + bpf[5] * -4.080258e-02
	    + bpf[6] * 3.353214e-02
	    + bpf[7] * -2.375732e-02
	    + bpf[8] * 8.203628e-03;

	/*
	 * The 1000/1200 sync signals are demodulated using a pair of
	 * quadrature multipliers and matched filters. However,
	 * synchronous demodulation at these frequencies is impractical,
	 * so only the signal amplitude is used. The I and Q quadrature
	 * sync signals are produced by multiplying the filtered signal
	 * by 1000-Hz (WWV) and 1200-Hz (WWVH) sine and cosine signals,
	 * respectively. The WWV and WWVH signals are processed by 800-
	 * ms synchronous matched filters and combined to produce the
	 * minute sync signal and detect which one (or both) the WWV or
	 * WWVH signal is present. The WWV and WWVH signals are also
	 * processed by 5-ms synchronous matched filters and combined to
	 * produce the second sync signal. The signals are scaled to
	 * produce unit energy at the maximum value.
	 *
	 * Note the master timing ramps, which run continuously. The
	 * minute counter (mphase) counts the samples in the minute,
	 * while the second counter (epoch) counts the samples in the
	 * second.
	 */
	up->mphase = (up->mphase + 1) % MINUTE;
	epoch = up->mphase % SECOND;

	/*
	 * WWV
	 */
	i = csinptr;
	csinptr = (csinptr + IN1000) % 80;

	dtemp = sintab[i] * syncx / (MS / 2.);
	ciamp -= cibuf[jptr];
	cibuf[jptr] = dtemp;
	ciamp += dtemp;
	csiamp -= csibuf[kptr];
	csibuf[kptr] = dtemp;
	csiamp += dtemp;

	i = (i + 20) % 80;
	dtemp = sintab[i] * syncx / (MS / 2.);
	cqamp -= cqbuf[jptr];
	cqbuf[jptr] = dtemp;
	cqamp += dtemp;
	csqamp -= csqbuf[kptr];
	csqbuf[kptr] = dtemp;
	csqamp += dtemp;

	sp = &up->mitig[up->achan].wwv;
	sp->amp = sqrt(ciamp * ciamp + cqamp * cqamp) / SYNCYC;
	if (!(up->status & MSYNC))
		wwv_qrz(peer, sp, (int)(pp->fudgetime1 * SECOND));

	/*
	 * WWVH
	 */
	i = hsinptr;
	hsinptr = (hsinptr + IN1200) % 80;

	dtemp = sintab[i] * syncx / (MS / 2.);
	hiamp -= hibuf[jptr];
	hibuf[jptr] = dtemp;
	hiamp += dtemp;
	hsiamp -= hsibuf[kptr];
	hsibuf[kptr] = dtemp;
	hsiamp += dtemp;

	i = (i + 20) % 80;
	dtemp = sintab[i] * syncx / (MS / 2.);
	hqamp -= hqbuf[jptr];
	hqbuf[jptr] = dtemp;
	hqamp += dtemp;
	hsqamp -= hsqbuf[kptr];
	hsqbuf[kptr] = dtemp;
	hsqamp += dtemp;

	rp = &up->mitig[up->achan].wwvh;
	rp->amp = sqrt(hiamp * hiamp + hqamp * hqamp) / SYNCYC;
	if (!(up->status & MSYNC))
		wwv_qrz(peer, rp, (int)(pp->fudgetime2 * SECOND));
	jptr = (jptr + 1) % SYNSIZ;
	kptr = (kptr + 1) % TCKSIZ;

	/*
	 * The following section is called once per minute. It does
	 * housekeeping and timeout functions and empties the dustbins.
	 */
	if (up->mphase == 0) {
		up->watch++;
		if (!(up->status & MSYNC)) {

			/*
			 * If minute sync has not been acquired before
			 * timeout, or if no signal is heard, the
			 * program cycles to the next frequency and
			 * tries again.
			 */
			if (!wwv_newchan(peer) || up->watch > ACQSN) {
#ifdef ICOM
				if (up->fd_icom > 0)
					up->dchan = (up->dchan + 1) %
					    NCHAN;
#endif /* ICOM */
				wwv_newgame(peer, up->dchan);
			}
		} else {

			/*
			 * If the leap bit is set, set the minute epoch
			 * back one second so the station processes
			 * don't miss a beat.
			 */
			if (up->status & LEPSEC) {
				up->mphase -= SECOND;
				if (up->mphase < 0)
					up->mphase += MINUTE;
			}
		}
	}

	/*
	 * When the channel metric reaches threshold and the second
	 * counter matches the minute epoch within the second, the
	 * driver has synchronized to the station. The second number is
	 * the remaining seconds until the next minute epoch, while the
	 * sync epoch is zero. Watch out for the first second; if
	 * already synchronized to the second, the buffered sync epoch
	 * must be set.
	 */
	if (up->status & MSYNC) {
		wwv_epoch(peer);
	} else if (up->sptr != NULL) {
		struct chan *cp;

		sp = up->sptr;
		if (sp->metric >= TTHR && epoch == sp->mepoch % SECOND)
		    {
			up->rsec = (60 - sp->mepoch / SECOND) % 60;
			up->rphase = 0;
			up->status |= MSYNC;
			up->watch = 0;
			if (!(up->status & SSYNC))
				up->repoch = up->yepoch = epoch;
			else
				up->repoch = up->yepoch;

			/*
			 * Clear the crud and initialize fairly.
			 */
			for (i = 0; i < NCHAN; i++) {
				cp = &up->mitig[i];
				cp->wwv.count = cp->wwv.reach = 0;
				cp->wwvh.count = cp->wwvh.reach = 0;
			}
			sp->count = sp->reach = 1;
		}
	}

	/*
	 * The second sync pulse is extracted using 5-ms (40 sample) FIR
	 * matched filters at 1000 Hz for WWV or 1200 Hz for WWVH. This
	 * pulse is used for the most precise synchronization, since if
	 * provides a resolution of one sample (125 us). The filters run
	 * only if the station has been reliably determined.
	 */
	if (up->status & SELV) {
		pdelay = (int)(pp->fudgetime1 * SECOND);
		mfsync = sqrt(csiamp * csiamp + csqamp * csqamp) /
		    TCKCYC;
	} else if (up->status & SELH) {
		pdelay = (int)(pp->fudgetime2 * SECOND);
		mfsync = sqrt(hsiamp * hsiamp + hsqamp * hsqamp) /
		    TCKCYC;
	} else {
		pdelay = 0;
		mfsync = 0;
	}

	/*
	 * Enhance the seconds sync pulse using a 1-s (8000-sample) comb
	 * filter. Correct for the FIR matched filter delay, which is 5
	 * ms for both the WWV and WWVH filters, and also for the
	 * propagation delay. Once each second look for second sync. If
	 * not in minute sync, fiddle the codec gain. Note the SNR is
	 * computed from the maximum sample and the envelope of the
	 * sample 10 ms before it, so if we slip more than a cycle the
	 * SNR should plummet. The signal is scaled to produce unit
	 * energy at the maximum value.
	 */
	dtemp = (epobuf[epoch] += (mfsync - epobuf[epoch]) /
	    up->avgint);
	if (dtemp > epomax) {
		epomax = dtemp;
		epopos = epoch;
	}
	if (epoch == 0) {
		int j;

		up->epomax = epomax;
		dtemp = 0;
		j = epopos - 10 * MS;
		if (j < 0)
			j += SECOND;
		up->eposnr = wwv_snr(epomax, epobuf[j]);
		epopos -= pdelay + TCKCYC * MS;
		if (epopos < 0)
			epopos += SECOND;
		wwv_endpoc(peer, epopos);
		if (!(up->status & SSYNC))
			up->alarm |= SYNERR;
		epomax = 0;
		if (!(up->status & MSYNC))
			wwv_gain(peer);
	}
}


/*
 * wwv_qrz - identify and acquire WWV/WWVH minute sync pulse