Home | History | Annotate | Line # | Download | only in ntpd
refclock_irig.c revision 1.5.6.2
      1 /*	$NetBSD: refclock_irig.c,v 1.5.6.2 2015/01/07 12:13:27 msaitoh Exp $	*/
      2 
      3 /*
      4  * refclock_irig - audio IRIG-B/E demodulator/decoder
      5  */
      6 #ifdef HAVE_CONFIG_H
      7 #include <config.h>
      8 #endif
      9 
     10 #if defined(REFCLOCK) && defined(CLOCK_IRIG)
     11 
     12 #include "ntpd.h"
     13 #include "ntp_io.h"
     14 #include "ntp_refclock.h"
     15 #include "ntp_calendar.h"
     16 #include "ntp_stdlib.h"
     17 
     18 #include <stdio.h>
     19 #include <ctype.h>
     20 #include <math.h>
     21 #ifdef HAVE_SYS_IOCTL_H
     22 #include <sys/ioctl.h>
     23 #endif /* HAVE_SYS_IOCTL_H */
     24 
     25 #include "audio.h"
     26 
     27 /*
     28  * Audio IRIG-B/E demodulator/decoder
     29  *
     30  * This driver synchronizes the computer time using data encoded in
     31  * IRIG-B/E signals commonly produced by GPS receivers and other timing
     32  * devices. The IRIG signal is an amplitude-modulated carrier with
     33  * pulse-width modulated data bits. For IRIG-B, the carrier frequency is
     34  * 1000 Hz and bit rate 100 b/s; for IRIG-E, the carrier frequenchy is
     35  * 100 Hz and bit rate 10 b/s. The driver automatically recognizes which
     36  & format is in use.
     37  *
     38  * The driver requires an audio codec or sound card with sampling rate 8
     39  * kHz and mu-law companding. This is the same standard as used by the
     40  * telephone industry and is supported by most hardware and operating
     41  * systems, including Solaris, SunOS, FreeBSD, NetBSD and Linux. In this
     42  * implementation, only one audio driver and codec can be supported on a
     43  * single machine.
     44  *
     45  * The program processes 8000-Hz mu-law companded samples using separate
     46  * signal filters for IRIG-B and IRIG-E, a comb filter, envelope
     47  * detector and automatic threshold corrector. Cycle crossings relative
     48  * to the corrected slice level determine the width of each pulse and
     49  * its value - zero, one or position identifier.
     50  *
     51  * The data encode 20 BCD digits which determine the second, minute,
     52  * hour and day of the year and sometimes the year and synchronization
     53  * condition. The comb filter exponentially averages the corresponding
     54  * samples of successive baud intervals in order to reliably identify
     55  * the reference carrier cycle. A type-II phase-lock loop (PLL) performs
     56  * additional integration and interpolation to accurately determine the
     57  * zero crossing of that cycle, which determines the reference
     58  * timestamp. A pulse-width discriminator demodulates the data pulses,
     59  * which are then encoded as the BCD digits of the timecode.
     60  *
     61  * The timecode and reference timestamp are updated once each second
     62  * with IRIG-B (ten seconds with IRIG-E) and local clock offset samples
     63  * saved for later processing. At poll intervals of 64 s, the saved
     64  * samples are processed by a trimmed-mean filter and used to update the
     65  * system clock.
     66  *
     67  * An automatic gain control feature provides protection against
     68  * overdriven or underdriven input signal amplitudes. It is designed to
     69  * maintain adequate demodulator signal amplitude while avoiding
     70  * occasional noise spikes. In order to assure reliable capture, the
     71  * decompanded input signal amplitude must be greater than 100 units and
     72  * the codec sample frequency error less than 250 PPM (.025 percent).
     73  *
     74  * Monitor Data
     75  *
     76  * The timecode format used for debugging and data recording includes
     77  * data helpful in diagnosing problems with the IRIG signal and codec
     78  * connections. The driver produces one line for each timecode in the
     79  * following format:
     80  *
     81  * 00 00 98 23 19:26:52 2782 143 0.694 10 0.3 66.5 3094572411.00027
     82  *
     83  * If clockstats is enabled, the most recent line is written to the
     84  * clockstats file every 64 s. If verbose recording is enabled (fudge
     85  * flag 4) each line is written as generated.
     86  *
     87  * The first field containes the error flags in hex, where the hex bits
     88  * are interpreted as below. This is followed by the year of century,
     89  * day of year and time of day. Note that the time of day is for the
     90  * previous minute, not the current time. The status indicator and year
     91  * are not produced by some IRIG devices and appear as zeros. Following
     92  * these fields are the carrier amplitude (0-3000), codec gain (0-255),
     93  * modulation index (0-1), time constant (4-10), carrier phase error
     94  * +-.5) and carrier frequency error (PPM). The last field is the on-
     95  * time timestamp in NTP format.
     96  *
     97  * The error flags are defined as follows in hex:
     98  *
     99  * x01	Low signal. The carrier amplitude is less than 100 units. This
    100  *	is usually the result of no signal or wrong input port.
    101  * x02	Frequency error. The codec frequency error is greater than 250
    102  *	PPM. This may be due to wrong signal format or (rarely)
    103  *	defective codec.
    104  * x04	Modulation error. The IRIG modulation index is less than 0.5.
    105  *	This is usually the result of an overdriven codec, wrong signal
    106  *	format or wrong input port.
    107  * x08	Frame synch error. The decoder frame does not match the IRIG
    108  *	frame. This is usually the result of an overdriven codec, wrong
    109  *	signal format or noisy IRIG signal. It may also be the result of
    110  *	an IRIG signature check which indicates a failure of the IRIG
    111  *	signal synchronization source.
    112  * x10	Data bit error. The data bit length is out of tolerance. This is
    113  *	usually the result of an overdriven codec, wrong signal format
    114  *	or noisy IRIG signal.
    115  * x20	Seconds numbering discrepancy. The decoder second does not match
    116  *	the IRIG second. This is usually the result of an overdriven
    117  *	codec, wrong signal format or noisy IRIG signal.
    118  * x40	Codec error (overrun). The machine is not fast enough to keep up
    119  *	with the codec.
    120  * x80	Device status error (Spectracom).
    121  *
    122  *
    123  * Once upon a time, an UltrSPARC 30 and Solaris 2.7 kept the clock
    124  * within a few tens of microseconds relative to the IRIG-B signal.
    125  * Accuracy with IRIG-E was about ten times worse. Unfortunately, Sun
    126  * broke the 2.7 audio driver in 2.8, which has a 10-ms sawtooth
    127  * modulation.
    128  *
    129  * Unlike other drivers, which can have multiple instantiations, this
    130  * one supports only one. It does not seem likely that more than one
    131  * audio codec would be useful in a single machine. More than one would
    132  * probably chew up too much CPU time anyway.
    133  *
    134  * Fudge factors
    135  *
    136  * Fudge flag4 causes the dubugging output described above to be
    137  * recorded in the clockstats file. Fudge flag2 selects the audio input
    138  * port, where 0 is the mike port (default) and 1 is the line-in port.
    139  * It does not seem useful to select the compact disc player port. Fudge
    140  * flag3 enables audio monitoring of the input signal. For this purpose,
    141  * the monitor gain is set t a default value. Fudgetime2 is used as a
    142  * frequency vernier for broken codec sample frequency.
    143  *
    144  * Alarm codes
    145  *
    146  * CEVNT_BADTIME	invalid date or time
    147  * CEVNT_TIMEOUT	no IRIG data since last poll
    148  */
    149 /*
    150  * Interface definitions
    151  */
    152 #define	DEVICE_AUDIO	"/dev/audio" /* audio device name */
    153 #define	PRECISION	(-17)	/* precision assumed (about 10 us) */
    154 #define	REFID		"IRIG"	/* reference ID */
    155 #define	DESCRIPTION	"Generic IRIG Audio Driver" /* WRU */
    156 #define	AUDIO_BUFSIZ	320	/* audio buffer size (40 ms) */
    157 #define SECOND		8000	/* nominal sample rate (Hz) */
    158 #define BAUD		80	/* samples per baud interval */
    159 #define OFFSET		128	/* companded sample offset */
    160 #define SIZE		256	/* decompanding table size */
    161 #define CYCLE		8	/* samples per bit */
    162 #define SUBFLD		10	/* bits per frame */
    163 #define FIELD		100	/* bits per second */
    164 #define MINTC		2	/* min PLL time constant */
    165 #define MAXTC		10	/* max PLL time constant max */
    166 #define	MAXAMP		3000.	/* maximum signal amplitude */
    167 #define	MINAMP		2000.	/* minimum signal amplitude */
    168 #define DRPOUT		100.	/* dropout signal amplitude */
    169 #define MODMIN		0.5	/* minimum modulation index */
    170 #define MAXFREQ		(250e-6 * SECOND) /* freq tolerance (.025%) */
    171 
    172 /*
    173  * The on-time synchronization point is the positive-going zero crossing
    174  * of the first cycle of the second. The IIR baseband filter phase delay
    175  * is 1.03 ms for IRIG-B and 3.47 ms for IRIG-E. The fudge value 2.68 ms
    176  * due to the codec and other causes was determined by calibrating to a
    177  * PPS signal from a GPS receiver.
    178  *
    179  * The results with a 2.4-GHz P4 running FreeBSD 6.1 are generally
    180  * within .02 ms short-term with .02 ms jitter. The processor load due
    181  * to the driver is 0.51 percent.
    182  */
    183 #define IRIG_B	((1.03 + 2.68) / 1000)	/* IRIG-B system delay (s) */
    184 #define IRIG_E	((3.47 + 2.68) / 1000)	/* IRIG-E system delay (s) */
    185 
    186 /*
    187  * Data bit definitions
    188  */
    189 #define BIT0		0	/* zero */
    190 #define BIT1		1	/* one */
    191 #define BITP		2	/* position identifier */
    192 
    193 /*
    194  * Error flags
    195  */
    196 #define IRIG_ERR_AMP	0x01	/* low carrier amplitude */
    197 #define IRIG_ERR_FREQ	0x02	/* frequency tolerance exceeded */
    198 #define IRIG_ERR_MOD	0x04	/* low modulation index */
    199 #define IRIG_ERR_SYNCH	0x08	/* frame synch error */
    200 #define IRIG_ERR_DECODE	0x10	/* frame decoding error */
    201 #define IRIG_ERR_CHECK	0x20	/* second numbering discrepancy */
    202 #define IRIG_ERR_ERROR	0x40	/* codec error (overrun) */
    203 #define IRIG_ERR_SIGERR	0x80	/* IRIG status error (Spectracom) */
    204 
    205 static	char	hexchar[] = "0123456789abcdef";
    206 
    207 /*
    208  * IRIG unit control structure
    209  */
    210 struct irigunit {
    211 	u_char	timecode[2 * SUBFLD + 1]; /* timecode string */
    212 	l_fp	timestamp;	/* audio sample timestamp */
    213 	l_fp	tick;		/* audio sample increment */
    214 	l_fp	refstamp;	/* reference timestamp */
    215 	l_fp	chrstamp;	/* baud timestamp */
    216 	l_fp	prvstamp;	/* previous baud timestamp */
    217 	double	integ[BAUD];	/* baud integrator */
    218 	double	phase, freq;	/* logical clock phase and frequency */
    219 	double	zxing;		/* phase detector integrator */
    220 	double	yxing;		/* cycle phase */
    221 	double	exing;		/* envelope phase */
    222 	double	modndx;		/* modulation index */
    223 	double	irig_b;		/* IRIG-B signal amplitude */
    224 	double	irig_e;		/* IRIG-E signal amplitude */
    225 	int	errflg;		/* error flags */
    226 	/*
    227 	 * Audio codec variables
    228 	 */
    229 	double	comp[SIZE];	/* decompanding table */
    230 	double	signal;		/* peak signal for AGC */
    231 	int	port;		/* codec port */
    232 	int	gain;		/* codec gain */
    233 	int	mongain;	/* codec monitor gain */
    234 	int	seccnt;		/* second interval counter */
    235 
    236 	/*
    237 	 * RF variables
    238 	 */
    239 	double	bpf[9];		/* IRIG-B filter shift register */
    240 	double	lpf[5];		/* IRIG-E filter shift register */
    241 	double	envmin, envmax;	/* envelope min and max */
    242 	double	slice;		/* envelope slice level */
    243 	double	intmin, intmax;	/* integrated envelope min and max */
    244 	double	maxsignal;	/* integrated peak amplitude */
    245 	double	noise;		/* integrated noise amplitude */
    246 	double	lastenv[CYCLE];	/* last cycle amplitudes */
    247 	double	lastint[CYCLE];	/* last integrated cycle amplitudes */
    248 	double	lastsig;	/* last carrier sample */
    249 	double	fdelay;		/* filter delay */
    250 	int	decim;		/* sample decimation factor */
    251 	int	envphase;	/* envelope phase */
    252 	int	envptr;		/* envelope phase pointer */
    253 	int	envsw;		/* envelope state */
    254 	int	envxing;	/* envelope slice crossing */
    255 	int	tc;		/* time constant */
    256 	int	tcount;		/* time constant counter */
    257 	int	badcnt;		/* decimation interval counter */
    258 
    259 	/*
    260 	 * Decoder variables
    261 	 */
    262 	int	pulse;		/* cycle counter */
    263 	int	cycles;		/* carrier cycles */
    264 	int	dcycles;	/* data cycles */
    265 	int	lastbit;	/* last code element */
    266 	int	second;		/* previous second */
    267 	int	bitcnt;		/* bit count in frame */
    268 	int	frmcnt;		/* bit count in second */
    269 	int	xptr;		/* timecode pointer */
    270 	int	bits;		/* demodulated bits */
    271 };
    272 
    273 /*
    274  * Function prototypes
    275  */
    276 static	int	irig_start	(int, struct peer *);
    277 static	void	irig_shutdown	(int, struct peer *);
    278 static	void	irig_receive	(struct recvbuf *);
    279 static	void	irig_poll	(int, struct peer *);
    280 
    281 /*
    282  * More function prototypes
    283  */
    284 static	void	irig_base	(struct peer *, double);
    285 static	void	irig_rf		(struct peer *, double);
    286 static	void	irig_baud	(struct peer *, int);
    287 static	void	irig_decode	(struct peer *, int);
    288 static	void	irig_gain	(struct peer *);
    289 
    290 /*
    291  * Transfer vector
    292  */
    293 struct	refclock refclock_irig = {
    294 	irig_start,		/* start up driver */
    295 	irig_shutdown,		/* shut down driver */
    296 	irig_poll,		/* transmit poll message */
    297 	noentry,		/* not used (old irig_control) */
    298 	noentry,		/* initialize driver (not used) */
    299 	noentry,		/* not used (old irig_buginfo) */
    300 	NOFLAGS			/* not used */
    301 };
    302 
    303 
    304 /*
    305  * irig_start - open the devices and initialize data for processing
    306  */
    307 static int
    308 irig_start(
    309 	int	unit,		/* instance number (used for PCM) */
    310 	struct peer *peer	/* peer structure pointer */
    311 	)
    312 {
    313 	struct refclockproc *pp;
    314 	struct irigunit *up;
    315 
    316 	/*
    317 	 * Local variables
    318 	 */
    319 	int	fd;		/* file descriptor */
    320 	int	i;		/* index */
    321 	double	step;		/* codec adjustment */
    322 
    323 	/*
    324 	 * Open audio device
    325 	 */
    326 	fd = audio_init(DEVICE_AUDIO, AUDIO_BUFSIZ, unit);
    327 	if (fd < 0)
    328 		return (0);
    329 #ifdef DEBUG
    330 	if (debug)
    331 		audio_show();
    332 #endif
    333 
    334 	/*
    335 	 * Allocate and initialize unit structure
    336 	 */
    337 	up = emalloc_zero(sizeof(*up));
    338 	pp = peer->procptr;
    339 	pp->io.clock_recv = irig_receive;
    340 	pp->io.srcclock = peer;
    341 	pp->io.datalen = 0;
    342 	pp->io.fd = fd;
    343 	if (!io_addclock(&pp->io)) {
    344 		close(fd);
    345 		pp->io.fd = -1;
    346 		free(up);
    347 		return (0);
    348 	}
    349 	pp->unitptr = up;
    350 
    351 	/*
    352 	 * Initialize miscellaneous variables
    353 	 */
    354 	peer->precision = PRECISION;
    355 	pp->clockdesc = DESCRIPTION;
    356 	memcpy((char *)&pp->refid, REFID, 4);
    357 	up->tc = MINTC;
    358 	up->decim = 1;
    359 	up->gain = 127;
    360 
    361 	/*
    362 	 * The companded samples are encoded sign-magnitude. The table
    363 	 * contains all the 256 values in the interest of speed.
    364 	 */
    365 	up->comp[0] = up->comp[OFFSET] = 0.;
    366 	up->comp[1] = 1; up->comp[OFFSET + 1] = -1.;
    367 	up->comp[2] = 3; up->comp[OFFSET + 2] = -3.;
    368 	step = 2.;
    369 	for (i = 3; i < OFFSET; i++) {
    370 		up->comp[i] = up->comp[i - 1] + step;
    371 		up->comp[OFFSET + i] = -up->comp[i];
    372 		if (i % 16 == 0)
    373 			step *= 2.;
    374 	}
    375 	DTOLFP(1. / SECOND, &up->tick);
    376 	return (1);
    377 }
    378 
    379 
    380 /*
    381  * irig_shutdown - shut down the clock
    382  */
    383 static void
    384 irig_shutdown(
    385 	int	unit,		/* instance number (not used) */
    386 	struct peer *peer	/* peer structure pointer */
    387 	)
    388 {
    389 	struct refclockproc *pp;
    390 	struct irigunit *up;
    391 
    392 	pp = peer->procptr;
    393 	up = pp->unitptr;
    394 	if (-1 != pp->io.fd)
    395 		io_closeclock(&pp->io);
    396 	if (NULL != up)
    397 		free(up);
    398 }
    399 
    400 
    401 /*
    402  * irig_receive - receive data from the audio device
    403  *
    404  * This routine reads input samples and adjusts the logical clock to
    405  * track the irig clock by dropping or duplicating codec samples.
    406  */
    407 static void
    408 irig_receive(
    409 	struct recvbuf *rbufp	/* receive buffer structure pointer */
    410 	)
    411 {
    412 	struct peer *peer;
    413 	struct refclockproc *pp;
    414 	struct irigunit *up;
    415 
    416 	/*
    417 	 * Local variables
    418 	 */
    419 	double	sample;		/* codec sample */
    420 	u_char	*dpt;		/* buffer pointer */
    421 	int	bufcnt;		/* buffer counter */
    422 	l_fp	ltemp;		/* l_fp temp */
    423 
    424 	peer = rbufp->recv_peer;
    425 	pp = peer->procptr;
    426 	up = pp->unitptr;
    427 
    428 	/*
    429 	 * Main loop - read until there ain't no more. Note codec
    430 	 * samples are bit-inverted.
    431 	 */
    432 	DTOLFP((double)rbufp->recv_length / SECOND, &ltemp);
    433 	L_SUB(&rbufp->recv_time, &ltemp);
    434 	up->timestamp = rbufp->recv_time;
    435 	dpt = rbufp->recv_buffer;
    436 	for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) {
    437 		sample = up->comp[~*dpt++ & 0xff];
    438 
    439 		/*
    440 		 * Variable frequency oscillator. The codec oscillator
    441 		 * runs at the nominal rate of 8000 samples per second,
    442 		 * or 125 us per sample. A frequency change of one unit
    443 		 * results in either duplicating or deleting one sample
    444 		 * per second, which results in a frequency change of
    445 		 * 125 PPM.
    446 		 */
    447 		up->phase += (up->freq + clock_codec) / SECOND;
    448 		up->phase += pp->fudgetime2 / 1e6;
    449 		if (up->phase >= .5) {
    450 			up->phase -= 1.;
    451 		} else if (up->phase < -.5) {
    452 			up->phase += 1.;
    453 			irig_rf(peer, sample);
    454 			irig_rf(peer, sample);
    455 		} else {
    456 			irig_rf(peer, sample);
    457 		}
    458 		L_ADD(&up->timestamp, &up->tick);
    459 		sample = fabs(sample);
    460 		if (sample > up->signal)
    461 			up->signal = sample;
    462 		up->signal += (sample - up->signal) /
    463 		    1000;
    464 
    465 		/*
    466 		 * Once each second, determine the IRIG format and gain.
    467 		 */
    468 		up->seccnt = (up->seccnt + 1) % SECOND;
    469 		if (up->seccnt == 0) {
    470 			if (up->irig_b > up->irig_e) {
    471 				up->decim = 1;
    472 				up->fdelay = IRIG_B;
    473 			} else {
    474 				up->decim = 10;
    475 				up->fdelay = IRIG_E;
    476 			}
    477 			up->irig_b = up->irig_e = 0;
    478 			irig_gain(peer);
    479 
    480 		}
    481 	}
    482 
    483 	/*
    484 	 * Set the input port and monitor gain for the next buffer.
    485 	 */
    486 	if (pp->sloppyclockflag & CLK_FLAG2)
    487 		up->port = 2;
    488 	else
    489 		up->port = 1;
    490 	if (pp->sloppyclockflag & CLK_FLAG3)
    491 		up->mongain = MONGAIN;
    492 	else
    493 		up->mongain = 0;
    494 }
    495 
    496 
    497 /*
    498  * irig_rf - RF processing
    499  *
    500  * This routine filters the RF signal using a bandass filter for IRIG-B
    501  * and a lowpass filter for IRIG-E. In case of IRIG-E, the samples are
    502  * decimated by a factor of ten. Note that the codec filters function as
    503  * roofing filters to attenuate both the high and low ends of the
    504  * passband. IIR filter coefficients were determined using Matlab Signal
    505  * Processing Toolkit.
    506  */
    507 static void
    508 irig_rf(
    509 	struct peer *peer,	/* peer structure pointer */
    510 	double	sample		/* current signal sample */
    511 	)
    512 {
    513 	struct refclockproc *pp;
    514 	struct irigunit *up;
    515 
    516 	/*
    517 	 * Local variables
    518 	 */
    519 	double	irig_b, irig_e;	/* irig filter outputs */
    520 
    521 	pp = peer->procptr;
    522 	up = pp->unitptr;
    523 
    524 	/*
    525 	 * IRIG-B filter. Matlab 4th-order IIR elliptic, 800-1200 Hz
    526 	 * bandpass, 0.3 dB passband ripple, -50 dB stopband ripple,
    527 	 * phase delay 1.03 ms.
    528 	 */
    529 	irig_b = (up->bpf[8] = up->bpf[7]) * 6.505491e-001;
    530 	irig_b += (up->bpf[7] = up->bpf[6]) * -3.875180e+000;
    531 	irig_b += (up->bpf[6] = up->bpf[5]) * 1.151180e+001;
    532 	irig_b += (up->bpf[5] = up->bpf[4]) * -2.141264e+001;
    533 	irig_b += (up->bpf[4] = up->bpf[3]) * 2.712837e+001;
    534 	irig_b += (up->bpf[3] = up->bpf[2]) * -2.384486e+001;
    535 	irig_b += (up->bpf[2] = up->bpf[1]) * 1.427663e+001;
    536 	irig_b += (up->bpf[1] = up->bpf[0]) * -5.352734e+000;
    537 	up->bpf[0] = sample - irig_b;
    538 	irig_b = up->bpf[0] * 4.952157e-003
    539 	    + up->bpf[1] * -2.055878e-002
    540 	    + up->bpf[2] * 4.401413e-002
    541 	    + up->bpf[3] * -6.558851e-002
    542 	    + up->bpf[4] * 7.462108e-002
    543 	    + up->bpf[5] * -6.558851e-002
    544 	    + up->bpf[6] * 4.401413e-002
    545 	    + up->bpf[7] * -2.055878e-002
    546 	    + up->bpf[8] * 4.952157e-003;
    547 	up->irig_b += irig_b * irig_b;
    548 
    549 	/*
    550 	 * IRIG-E filter. Matlab 4th-order IIR elliptic, 130-Hz lowpass,
    551 	 * 0.3 dB passband ripple, -50 dB stopband ripple, phase delay
    552 	 * 3.47 ms.
    553 	 */
    554 	irig_e = (up->lpf[4] = up->lpf[3]) * 8.694604e-001;
    555 	irig_e += (up->lpf[3] = up->lpf[2]) * -3.589893e+000;
    556 	irig_e += (up->lpf[2] = up->lpf[1]) * 5.570154e+000;
    557 	irig_e += (up->lpf[1] = up->lpf[0]) * -3.849667e+000;
    558 	up->lpf[0] = sample - irig_e;
    559 	irig_e = up->lpf[0] * 3.215696e-003
    560 	    + up->lpf[1] * -1.174951e-002
    561 	    + up->lpf[2] * 1.712074e-002
    562 	    + up->lpf[3] * -1.174951e-002
    563 	    + up->lpf[4] * 3.215696e-003;
    564 	up->irig_e += irig_e * irig_e;
    565 
    566 	/*
    567 	 * Decimate by a factor of either 1 (IRIG-B) or 10 (IRIG-E).
    568 	 */
    569 	up->badcnt = (up->badcnt + 1) % up->decim;
    570 	if (up->badcnt == 0) {
    571 		if (up->decim == 1)
    572 			irig_base(peer, irig_b);
    573 		else
    574 			irig_base(peer, irig_e);
    575 	}
    576 }
    577 
    578 /*
    579  * irig_base - baseband processing
    580  *
    581  * This routine processes the baseband signal and demodulates the AM
    582  * carrier using a synchronous detector. It then synchronizes to the
    583  * data frame at the baud rate and decodes the width-modulated data
    584  * pulses.
    585  */
    586 static void
    587 irig_base(
    588 	struct peer *peer,	/* peer structure pointer */
    589 	double	sample		/* current signal sample */
    590 	)
    591 {
    592 	struct refclockproc *pp;
    593 	struct irigunit *up;
    594 
    595 	/*
    596 	 * Local variables
    597 	 */
    598 	double	lope;		/* integrator output */
    599 	double	env;		/* envelope detector output */
    600 	double	dtemp;
    601 	int	carphase;	/* carrier phase */
    602 
    603 	pp = peer->procptr;
    604 	up = pp->unitptr;
    605 
    606 	/*
    607 	 * Synchronous baud integrator. Corresponding samples of current
    608 	 * and past baud intervals are integrated to refine the envelope
    609 	 * amplitude and phase estimate. We keep one cycle (1 ms) of the
    610 	 * raw data and one baud (10 ms) of the integrated data.
    611 	 */
    612 	up->envphase = (up->envphase + 1) % BAUD;
    613 	up->integ[up->envphase] += (sample - up->integ[up->envphase]) /
    614 	    (5 * up->tc);
    615 	lope = up->integ[up->envphase];
    616 	carphase = up->envphase % CYCLE;
    617 	up->lastenv[carphase] = sample;
    618 	up->lastint[carphase] = lope;
    619 
    620 	/*
    621 	 * Phase detector. Find the negative-going zero crossing
    622 	 * relative to sample 4 in the 8-sample sycle. A phase change of
    623 	 * 360 degrees produces an output change of one unit.
    624 	 */
    625 	if (up->lastsig > 0 && lope <= 0)
    626 		up->zxing += (double)(carphase - 4) / CYCLE;
    627 	up->lastsig = lope;
    628 
    629 	/*
    630 	 * End of the baud. Update signal/noise estimates and PLL
    631 	 * phase, frequency and time constant.
    632 	 */
    633 	if (up->envphase == 0) {
    634 		up->maxsignal = up->intmax; up->noise = up->intmin;
    635 		up->intmin = 1e6; up->intmax = -1e6;
    636 		if (up->maxsignal < DRPOUT)
    637 			up->errflg |= IRIG_ERR_AMP;
    638 		if (up->maxsignal > 0)
    639 			up->modndx = (up->maxsignal - up->noise) /
    640 			    up->maxsignal;
    641  		else
    642 			up->modndx = 0;
    643 		if (up->modndx < MODMIN)
    644 			up->errflg |= IRIG_ERR_MOD;
    645 		if (up->errflg & (IRIG_ERR_AMP | IRIG_ERR_FREQ |
    646 		   IRIG_ERR_MOD | IRIG_ERR_SYNCH)) {
    647 			up->tc = MINTC;
    648 			up->tcount = 0;
    649 		}
    650 
    651 		/*
    652 		 * Update PLL phase and frequency. The PLL time constant
    653 		 * is set initially to stabilize the frequency within a
    654 		 * minute or two, then increases to the maximum. The
    655 		 * frequency is clamped so that the PLL capture range
    656 		 * cannot be exceeded.
    657 		 */
    658 		dtemp = up->zxing * up->decim / BAUD;
    659 		up->yxing = dtemp;
    660 		up->zxing = 0.;
    661 		up->phase += dtemp / up->tc;
    662 		up->freq += dtemp / (4. * up->tc * up->tc);
    663 		if (up->freq > MAXFREQ) {
    664 			up->freq = MAXFREQ;
    665 			up->errflg |= IRIG_ERR_FREQ;
    666 		} else if (up->freq < -MAXFREQ) {
    667 			up->freq = -MAXFREQ;
    668 			up->errflg |= IRIG_ERR_FREQ;
    669 		}
    670 	}
    671 
    672 	/*
    673 	 * Synchronous demodulator. There are eight samples in the cycle
    674 	 * and ten cycles in the baud. Since the PLL has aligned the
    675 	 * negative-going zero crossing at sample 4, the maximum
    676 	 * amplitude is at sample 2 and minimum at sample 6. The
    677 	 * beginning of the data pulse is determined from the integrated
    678 	 * samples, while the end of the pulse is determined from the
    679 	 * raw samples. The raw data bits are demodulated relative to
    680 	 * the slice level and left-shifted in the decoding register.
    681 	 */
    682 	if (carphase != 7)
    683 		return;
    684 
    685 	lope = (up->lastint[2] - up->lastint[6]) / 2.;
    686 	if (lope > up->intmax)
    687 		up->intmax = lope;
    688 	if (lope < up->intmin)
    689 		up->intmin = lope;
    690 
    691 	/*
    692 	 * Pulse code demodulator and reference timestamp. The decoder
    693 	 * looks for a sequence of ten bits; the first two bits must be
    694 	 * one, the last two bits must be zero. Frame synch is asserted
    695 	 * when three correct frames have been found.
    696 	 */
    697 	up->pulse = (up->pulse + 1) % 10;
    698 	up->cycles <<= 1;
    699 	if (lope >= (up->maxsignal + up->noise) / 2.)
    700 		up->cycles |= 1;
    701 	if ((up->cycles & 0x303c0f03) == 0x300c0300) {
    702 		if (up->pulse != 0)
    703 			up->errflg |= IRIG_ERR_SYNCH;
    704 		up->pulse = 0;
    705 	}
    706 
    707 	/*
    708 	 * Assemble the baud and max/min to get the slice level for the
    709 	 * next baud. The slice level is based on the maximum over the
    710 	 * first two bits and the minimum over the last two bits, with
    711 	 * the slice level halfway between the maximum and minimum.
    712 	 */
    713 	env = (up->lastenv[2] - up->lastenv[6]) / 2.;
    714 	up->dcycles <<= 1;
    715 	if (env >= up->slice)
    716 		up->dcycles |= 1;
    717 	switch(up->pulse) {
    718 
    719 	case 0:
    720 		irig_baud(peer, up->dcycles);
    721 		if (env < up->envmin)
    722 			up->envmin = env;
    723 		up->slice = (up->envmax + up->envmin) / 2;
    724 		up->envmin = 1e6; up->envmax = -1e6;
    725 		break;
    726 
    727 	case 1:
    728 		up->envmax = env;
    729 		break;
    730 
    731 	case 2:
    732 		if (env > up->envmax)
    733 			up->envmax = env;
    734 		break;
    735 
    736 	case 9:
    737 		up->envmin = env;
    738 		break;
    739 	}
    740 }
    741 
    742 /*
    743  * irig_baud - update the PLL and decode the pulse-width signal
    744  */
    745 static void
    746 irig_baud(
    747 	struct peer *peer,	/* peer structure pointer */
    748 	int	bits		/* decoded bits */
    749 	)
    750 {
    751 	struct refclockproc *pp;
    752 	struct irigunit *up;
    753 	double	dtemp;
    754 	l_fp	ltemp;
    755 
    756         pp = peer->procptr;
    757 	up = pp->unitptr;
    758 
    759 	/*
    760 	 * The PLL time constant starts out small, in order to
    761 	 * sustain a frequency tolerance of 250 PPM. It
    762 	 * gradually increases as the loop settles down. Note
    763 	 * that small wiggles are not believed, unless they
    764 	 * persist for lots of samples.
    765 	 */
    766 	up->exing = -up->yxing;
    767 	if (abs(up->envxing - up->envphase) <= 1) {
    768 		up->tcount++;
    769 		if (up->tcount > 20 * up->tc) {
    770 			up->tc++;
    771 			if (up->tc > MAXTC)
    772 				up->tc = MAXTC;
    773 			up->tcount = 0;
    774 			up->envxing = up->envphase;
    775 		} else {
    776 			up->exing -= up->envxing - up->envphase;
    777 		}
    778 	} else {
    779 		up->tcount = 0;
    780 		up->envxing = up->envphase;
    781 	}
    782 
    783 	/*
    784 	 * Strike the baud timestamp as the positive zero crossing of
    785 	 * the first bit, accounting for the codec delay and filter
    786 	 * delay.
    787 	 */
    788 	up->prvstamp = up->chrstamp;
    789 	dtemp = up->decim * (up->exing / SECOND) + up->fdelay;
    790 	DTOLFP(dtemp, &ltemp);
    791 	up->chrstamp = up->timestamp;
    792 	L_SUB(&up->chrstamp, &ltemp);
    793 
    794 	/*
    795 	 * The data bits are collected in ten-bit bauds. The first two
    796 	 * bits are not used. The resulting patterns represent runs of
    797 	 * 0-1 bits (0), 2-4 bits (1) and 5-7 bits (PI). The remaining
    798 	 * 8-bit run represents a soft error and is treated as 0.
    799 	 */
    800 	switch (up->dcycles & 0xff) {
    801 
    802 	case 0x00:		/* 0-1 bits (0) */
    803 	case 0x80:
    804 		irig_decode(peer, BIT0);
    805 		break;
    806 
    807 	case 0xc0:		/* 2-4 bits (1) */
    808 	case 0xe0:
    809 	case 0xf0:
    810 		irig_decode(peer, BIT1);
    811 		break;
    812 
    813 	case 0xf8:		/* (5-7 bits (PI) */
    814 	case 0xfc:
    815 	case 0xfe:
    816 		irig_decode(peer, BITP);
    817 		break;
    818 
    819 	default:		/* 8 bits (error) */
    820 		irig_decode(peer, BIT0);
    821 		up->errflg |= IRIG_ERR_DECODE;
    822 	}
    823 }
    824 
    825 
    826 /*
    827  * irig_decode - decode the data
    828  *
    829  * This routine assembles bauds into digits, digits into frames and
    830  * frames into the timecode fields. Bits can have values of zero, one
    831  * or position identifier. There are four bits per digit, ten digits per
    832  * frame and ten frames per second.
    833  */
    834 static void
    835 irig_decode(
    836 	struct	peer *peer,	/* peer structure pointer */
    837 	int	bit		/* data bit (0, 1 or 2) */
    838 	)
    839 {
    840 	struct refclockproc *pp;
    841 	struct irigunit *up;
    842 
    843 	/*
    844 	 * Local variables
    845 	 */
    846 	int	syncdig;	/* sync digit (Spectracom) */
    847 	char	sbs[6 + 1];	/* binary seconds since 0h */
    848 	char	spare[2 + 1];	/* mulligan digits */
    849 	int	temp;
    850 
    851 	syncdig = 0;
    852 	pp = peer->procptr;
    853 	up = pp->unitptr;
    854 
    855 	/*
    856 	 * Assemble frame bits.
    857 	 */
    858 	up->bits >>= 1;
    859 	if (bit == BIT1) {
    860 		up->bits |= 0x200;
    861 	} else if (bit == BITP && up->lastbit == BITP) {
    862 
    863 		/*
    864 		 * Frame sync - two adjacent position identifiers, which
    865 		 * mark the beginning of the second. The reference time
    866 		 * is the beginning of the second position identifier,
    867 		 * so copy the character timestamp to the reference
    868 		 * timestamp.
    869 		 */
    870 		if (up->frmcnt != 1)
    871 			up->errflg |= IRIG_ERR_SYNCH;
    872 		up->frmcnt = 1;
    873 		up->refstamp = up->prvstamp;
    874 	}
    875 	up->lastbit = bit;
    876 	if (up->frmcnt % SUBFLD == 0) {
    877 
    878 		/*
    879 		 * End of frame. Encode two hexadecimal digits in
    880 		 * little-endian timecode field. Note frame 1 is shifted
    881 		 * right one bit to account for the marker PI.
    882 		 */
    883 		temp = up->bits;
    884 		if (up->frmcnt == 10)
    885 			temp >>= 1;
    886 		if (up->xptr >= 2) {
    887 			up->timecode[--up->xptr] = hexchar[temp & 0xf];
    888 			up->timecode[--up->xptr] = hexchar[(temp >> 5) &
    889 			    0xf];
    890 		}
    891 		if (up->frmcnt == 0) {
    892 
    893 			/*
    894 			 * End of second. Decode the timecode and wind
    895 			 * the clock. Not all IRIG generators have the
    896 			 * year; if so, it is nonzero after year 2000.
    897 			 * Not all have the hardware status bit; if so,
    898 			 * it is lit when the source is okay and dim
    899 			 * when bad. We watch this only if the year is
    900 			 * nonzero. Not all are configured for signature
    901 			 * control. If so, all BCD digits are set to
    902 			 * zero if the source is bad. In this case the
    903 			 * refclock_process() will reject the timecode
    904 			 * as invalid.
    905 			 */
    906 			up->xptr = 2 * SUBFLD;
    907 			if (sscanf((char *)up->timecode,
    908 			   "%6s%2d%1d%2s%3d%2d%2d%2d", sbs, &pp->year,
    909 			    &syncdig, spare, &pp->day, &pp->hour,
    910 			    &pp->minute, &pp->second) != 8)
    911 				pp->leap = LEAP_NOTINSYNC;
    912 			else
    913 				pp->leap = LEAP_NOWARNING;
    914 			up->second = (up->second + up->decim) % 60;
    915 
    916 			/*
    917 			 * Raise an alarm if the day field is zero,
    918 			 * which happens when signature control is
    919 			 * enabled and the device has lost
    920 			 * synchronization. Raise an alarm if the year
    921 			 * field is nonzero and the sync indicator is
    922 			 * zero, which happens when a Spectracom radio
    923 			 * has lost synchronization. Raise an alarm if
    924 			 * the expected second does not agree with the
    925 			 * decoded second, which happens with a garbled
    926 			 * IRIG signal. We are very particular.
    927 			 */
    928 			if (pp->day == 0 || (pp->year != 0 && syncdig ==
    929 			    0))
    930 				up->errflg |= IRIG_ERR_SIGERR;
    931 			if (pp->second != up->second)
    932 				up->errflg |= IRIG_ERR_CHECK;
    933 			up->second = pp->second;
    934 
    935 			/*
    936 			 * Wind the clock only if there are no errors
    937 			 * and the time constant has reached the
    938 			 * maximum.
    939 			 */
    940 			if (up->errflg == 0 && up->tc == MAXTC) {
    941 				pp->lastref = pp->lastrec;
    942 				pp->lastrec = up->refstamp;
    943 				if (!refclock_process(pp))
    944 					refclock_report(peer,
    945 					    CEVNT_BADTIME);
    946 			}
    947 			snprintf(pp->a_lastcode, sizeof(pp->a_lastcode),
    948 			    "%02x %02d %03d %02d:%02d:%02d %4.0f %3d %6.3f %2d %6.2f %6.1f %s",
    949 			    up->errflg, pp->year, pp->day,
    950 			    pp->hour, pp->minute, pp->second,
    951 			    up->maxsignal, up->gain, up->modndx,
    952 			    up->tc, up->exing * 1e6 / SECOND, up->freq *
    953 			    1e6 / SECOND, ulfptoa(&pp->lastrec, 6));
    954 			pp->lencode = strlen(pp->a_lastcode);
    955 			up->errflg = 0;
    956 			if (pp->sloppyclockflag & CLK_FLAG4) {
    957 				record_clock_stats(&peer->srcadr,
    958 				    pp->a_lastcode);
    959 #ifdef DEBUG
    960 				if (debug)
    961 					printf("irig %s\n",
    962 					    pp->a_lastcode);
    963 #endif /* DEBUG */
    964 			}
    965 		}
    966 	}
    967 	up->frmcnt = (up->frmcnt + 1) % FIELD;
    968 }
    969 
    970 
    971 /*
    972  * irig_poll - called by the transmit procedure
    973  *
    974  * This routine sweeps up the timecode updates since the last poll. For
    975  * IRIG-B there should be at least 60 updates; for IRIG-E there should
    976  * be at least 6. If nothing is heard, a timeout event is declared.
    977  */
    978 static void
    979 irig_poll(
    980 	int	unit,		/* instance number (not used) */
    981 	struct peer *peer	/* peer structure pointer */
    982 	)
    983 {
    984 	struct refclockproc *pp;
    985 
    986 	pp = peer->procptr;
    987 
    988 	if (pp->coderecv == pp->codeproc) {
    989 		refclock_report(peer, CEVNT_TIMEOUT);
    990 		return;
    991 
    992 	}
    993 	refclock_receive(peer);
    994 	if (!(pp->sloppyclockflag & CLK_FLAG4)) {
    995 		record_clock_stats(&peer->srcadr, pp->a_lastcode);
    996 #ifdef DEBUG
    997 		if (debug)
    998 			printf("irig %s\n", pp->a_lastcode);
    999 #endif /* DEBUG */
   1000 	}
   1001 	pp->polls++;
   1002 
   1003 }
   1004 
   1005 
   1006 /*
   1007  * irig_gain - adjust codec gain
   1008  *
   1009  * This routine is called at the end of each second. It uses the AGC to
   1010  * bradket the maximum signal level between MINAMP and MAXAMP to avoid
   1011  * hunting. The routine also jiggles the input port and selectively
   1012  * mutes the monitor.
   1013  */
   1014 static void
   1015 irig_gain(
   1016 	struct peer *peer	/* peer structure pointer */
   1017 	)
   1018 {
   1019 	struct refclockproc *pp;
   1020 	struct irigunit *up;
   1021 
   1022 	pp = peer->procptr;
   1023 	up = pp->unitptr;
   1024 
   1025 	/*
   1026 	 * Apparently, the codec uses only the high order bits of the
   1027 	 * gain control field. Thus, it may take awhile for changes to
   1028 	 * wiggle the hardware bits.
   1029 	 */
   1030 	if (up->maxsignal < MINAMP) {
   1031 		up->gain += 4;
   1032 		if (up->gain > MAXGAIN)
   1033 			up->gain = MAXGAIN;
   1034 	} else if (up->maxsignal > MAXAMP) {
   1035 		up->gain -= 4;
   1036 		if (up->gain < 0)
   1037 			up->gain = 0;
   1038 	}
   1039 	audio_gain(up->gain, up->mongain, up->port);
   1040 }
   1041 
   1042 
   1043 #else
   1044 int refclock_irig_bs;
   1045 #endif /* REFCLOCK */
   1046