2  * refclock_irig - audio IRIG-B/E demodulator/decoder
26  * Audio IRIG-B/E demodulator/decoder
29  * IRIG-B/E signals commonly produced by GPS receivers and other timing
30  * devices. The IRIG signal is an amplitude-modulated carrier with
31  * pulse-width modulated data bits. For IRIG-B, the carrier frequency is
32  * 1000 Hz and bit rate 100 b/s; for IRIG-E, the carrier frequenchy is
37  * kHz and mu-law companding. This is the same standard as used by the
43  * The program processes 8000-Hz mu-law companded samples using separate
44  * signal filters for IRIG-B and IRIG-E, a comb filter, envelope
46  * to the corrected slice level determine the width of each pulse and
47  * its value - zero, one or position identifier.
53  * the reference carrier cycle. A type-II phase-lock loop (PLL) performs
56  * timestamp. A pulse-width discriminator demodulates the data pulses,
60  * with IRIG-B (ten seconds with IRIG-E) and local clock offset samples
62  * samples are processed by a trimmed-mean filter and used to update the
76  * connections. The driver produces one line for each timecode in the
81  * If clockstats is enabled, the most recent line is written to the
83  * flag 4) each line is written as generated.
90  * these fields are the carrier amplitude (0-3000), codec gain (0-255),
91  * modulation index (0-1), time constant (4-10), carrier phase error
92  * +-.5) and carrier frequency error (PPM). The last field is the on-
122  * within a few tens of microseconds relative to the IRIG-B signal.
123  * Accuracy with IRIG-E was about ten times worse. Unfortunately, Sun
124  * broke the 2.7 audio driver in 2.8, which has a 10-ms sawtooth
136  * port, where 0 is the mike port (default) and 1 is the line-in port.
151 #define	PRECISION	(-17)	/* precision assumed (about 10 us) */
156 #define BAUD		80	/* samples per baud interval */
159 #define CYCLE		8	/* samples per bit */
160 #define SUBFLD		10	/* bits per frame */
161 #define FIELD		100	/* bits per second */
168 #define MAXFREQ		(250e-6 * SECOND) /* freq tolerance (.025%) */
171  * The on-time synchronization point is the positive-going zero crossing
173  * is 1.03 ms for IRIG-B and 3.47 ms for IRIG-E. The fudge value 2.68 ms
177  * The results with a 2.4-GHz P4 running FreeBSD 6.1 are generally
178  * within .02 ms short-term with .02 ms jitter. The processor load due
181 #define IRIG_B	((1.03 + 2.68) / 1000)	/* IRIG-B system delay (s) */
182 #define IRIG_E	((3.47 + 2.68) / 1000)	/* IRIG-E system delay (s) */
221 	double	irig_b;		/* IRIG-B signal amplitude */
222 	double	irig_e;		/* IRIG-E signal amplitude */
237 	double	bpf[9];		/* IRIG-B filter shift register */
238 	double	lpf[5];		/* IRIG-E filter shift register */
240 	double	slice;		/* envelope slice level */  member
252 	int	envxing;	/* envelope slice crossing */
303  * irig_start - open the devices and initialize data for processing
336 	pp = peer->procptr;  in irig_start()
337 	pp->io.clock_recv = irig_receive;  in irig_start()
338 	pp->io.srcclock = peer;  in irig_start()
339 	pp->io.datalen = 0;  in irig_start()
340 	pp->io.fd = fd;  in irig_start()
341 	if (!io_addclock(&pp->io)) {  in irig_start()
343 		pp->io.fd = -1;  in irig_start()
347 	pp->unitptr = up;  in irig_start()
352 	peer->precision = PRECISION;  in irig_start()
353 	pp->clockdesc = DESCRIPTION;  in irig_start()
354 	memcpy((char *)&pp->refid, REFID, 4);  in irig_start()
355 	up->tc = MINTC;  in irig_start()
356 	up->decim = 1;  in irig_start()
357 	up->gain = 127;  in irig_start()
360 	 * The companded samples are encoded sign-magnitude. The table  in irig_start()
363 	up->comp[0] = up->comp[OFFSET] = 0.;  in irig_start()
364 	up->comp[1] = 1; up->comp[OFFSET + 1] = -1.;  in irig_start()
365 	up->comp[2] = 3; up->comp[OFFSET + 2] = -3.;  in irig_start()
368 		up->comp[i] = up->comp[i - 1] + step;  in irig_start()
369 		up->comp[OFFSET + i] = -up->comp[i];  in irig_start()
373 	DTOLFP(1. / SECOND, &up->tick);  in irig_start()
379  * irig_shutdown - shut down the clock
390 	pp = peer->procptr;  in irig_shutdown()
391 	up = pp->unitptr;  in irig_shutdown()
392 	if (-1 != pp->io.fd)  in irig_shutdown()
393 		io_closeclock(&pp->io);  in irig_shutdown()
400  * irig_receive - receive data from the audio device
422 	peer = rbufp->recv_peer;  in irig_receive()
423 	pp = peer->procptr;  in irig_receive()
424 	up = pp->unitptr;  in irig_receive()
427 	 * Main loop - read until there ain't no more. Note codec  in irig_receive()
428 	 * samples are bit-inverted.  in irig_receive()
430 	DTOLFP((double)rbufp->recv_length / SECOND, &ltemp);  in irig_receive()
431 	L_SUB(&rbufp->recv_time, &ltemp);  in irig_receive()
432 	up->timestamp = rbufp->recv_time;  in irig_receive()
433 	dpt = rbufp->recv_buffer;  in irig_receive()
434 	for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) {  in irig_receive()
435 		sample = up->comp[~*dpt++ & 0xff];  in irig_receive()
439 		 * runs at the nominal rate of 8000 samples per second,  in irig_receive()
440 		 * or 125 us per sample. A frequency change of one unit  in irig_receive()
442 		 * per second, which results in a frequency change of  in irig_receive()
445 		up->phase += (up->freq + clock_codec) / SECOND;  in irig_receive()
446 		up->phase += pp->fudgetime2 / 1e6;  in irig_receive()
447 		if (up->phase >= .5) {  in irig_receive()
448 			up->phase -= 1.;  in irig_receive()
449 		} else if (up->phase < -.5) {  in irig_receive()
450 			up->phase += 1.;  in irig_receive()
456 		L_ADD(&up->timestamp, &up->tick);  in irig_receive()
458 		if (sample > up->signal)  in irig_receive()
459 			up->signal = sample;  in irig_receive()
460 		up->signal += (sample - up->signal) /  in irig_receive()
466 		up->seccnt = (up->seccnt + 1) % SECOND;  in irig_receive()
467 		if (up->seccnt == 0) {  in irig_receive()
468 			if (up->irig_b > up->irig_e) {  in irig_receive()
469 				up->decim = 1;  in irig_receive()
470 				up->fdelay = IRIG_B;  in irig_receive()
472 				up->decim = 10;  in irig_receive()
473 				up->fdelay = IRIG_E;  in irig_receive()
475 			up->irig_b = up->irig_e = 0;  in irig_receive()
484 	if (pp->sloppyclockflag & CLK_FLAG2)  in irig_receive()
485 		up->port = 2;  in irig_receive()
487 		up->port = 1;  in irig_receive()
488 	if (pp->sloppyclockflag & CLK_FLAG3)  in irig_receive()
489 		up->mongain = MONGAIN;  in irig_receive()
491 		up->mongain = 0;  in irig_receive()
496  * irig_rf - RF processing
498  * This routine filters the RF signal using a bandass filter for IRIG-B
499  * and a lowpass filter for IRIG-E. In case of IRIG-E, the samples are
519 	pp = peer->procptr;  in irig_rf()
520 	up = pp->unitptr;  in irig_rf()
523 	 * IRIG-B filter. Matlab 4th-order IIR elliptic, 800-1200 Hz  in irig_rf()
524 	 * bandpass, 0.3 dB passband ripple, -50 dB stopband ripple,  in irig_rf()
527 	irig_b = (up->bpf[8] = up->bpf[7]) * 6.505491e-001;  in irig_rf()
528 	irig_b += (up->bpf[7] = up->bpf[6]) * -3.875180e+000;  in irig_rf()
529 	irig_b += (up->bpf[6] = up->bpf[5]) * 1.151180e+001;  in irig_rf()
530 	irig_b += (up->bpf[5] = up->bpf[4]) * -2.141264e+001;  in irig_rf()
531 	irig_b += (up->bpf[4] = up->bpf[3]) * 2.712837e+001;  in irig_rf()
532 	irig_b += (up->bpf[3] = up->bpf[2]) * -2.384486e+001;  in irig_rf()
533 	irig_b += (up->bpf[2] = up->bpf[1]) * 1.427663e+001;  in irig_rf()
534 	irig_b += (up->bpf[1] = up->bpf[0]) * -5.352734e+000;  in irig_rf()
535 	up->bpf[0] = sample - irig_b;  in irig_rf()
536 	irig_b = up->bpf[0] * 4.952157e-003  in irig_rf()
537 	    + up->bpf[1] * -2.055878e-002  in irig_rf()
538 	    + up->bpf[2] * 4.401413e-002  in irig_rf()
539 	    + up->bpf[3] * -6.558851e-002  in irig_rf()
540 	    + up->bpf[4] * 7.462108e-002  in irig_rf()
541 	    + up->bpf[5] * -6.558851e-002  in irig_rf()
542 	    + up->bpf[6] * 4.401413e-002  in irig_rf()
543 	    + up->bpf[7] * -2.055878e-002  in irig_rf()
544 	    + up->bpf[8] * 4.952157e-003;  in irig_rf()
545 	up->irig_b += irig_b * irig_b;  in irig_rf()
548 	 * IRIG-E filter. Matlab 4th-order IIR elliptic, 130-Hz lowpass,  in irig_rf()
549 	 * 0.3 dB passband ripple, -50 dB stopband ripple, phase delay  in irig_rf()
552 	irig_e = (up->lpf[4] = up->lpf[3]) * 8.694604e-001;  in irig_rf()
553 	irig_e += (up->lpf[3] = up->lpf[2]) * -3.589893e+000;  in irig_rf()
554 	irig_e += (up->lpf[2] = up->lpf[1]) * 5.570154e+000;  in irig_rf()
555 	irig_e += (up->lpf[1] = up->lpf[0]) * -3.849667e+000;  in irig_rf()
556 	up->lpf[0] = sample - irig_e;  in irig_rf()
557 	irig_e = up->lpf[0] * 3.215696e-003  in irig_rf()
558 	    + up->lpf[1] * -1.174951e-002  in irig_rf()
559 	    + up->lpf[2] * 1.712074e-002  in irig_rf()
560 	    + up->lpf[3] * -1.174951e-002  in irig_rf()
561 	    + up->lpf[4] * 3.215696e-003;  in irig_rf()
562 	up->irig_e += irig_e * irig_e;  in irig_rf()
565 	 * Decimate by a factor of either 1 (IRIG-B) or 10 (IRIG-E).  in irig_rf()
567 	up->badcnt = (up->badcnt + 1) % up->decim;  in irig_rf()
568 	if (up->badcnt == 0) {  in irig_rf()
569 		if (up->decim == 1)  in irig_rf()
577  * irig_base - baseband processing
581  * data frame at the baud rate and decodes the width-modulated data
601 	pp = peer->procptr;  in irig_base()
602 	up = pp->unitptr;  in irig_base()
610 	up->envphase = (up->envphase + 1) % BAUD;  in irig_base()
611 	up->integ[up->envphase] += (sample - up->integ[up->envphase]) /  in irig_base()
612 	    (5 * up->tc);  in irig_base()
613 	lope = up->integ[up->envphase];  in irig_base()
614 	carphase = up->envphase % CYCLE;  in irig_base()
615 	up->lastenv[carphase] = sample;  in irig_base()
616 	up->lastint[carphase] = lope;  in irig_base()
619 	 * Phase detector. Find the negative-going zero crossing  in irig_base()
620 	 * relative to sample 4 in the 8-sample sycle. A phase change of  in irig_base()
623 	if (up->lastsig > 0 && lope <= 0)  in irig_base()
624 		up->zxing += (double)(carphase - 4) / CYCLE;  in irig_base()
625 	up->lastsig = lope;  in irig_base()
631 	if (up->envphase == 0) {  in irig_base()
632 		up->maxsignal = up->intmax; up->noise = up->intmin;  in irig_base()
633 		up->intmin = 1e6; up->intmax = -1e6;  in irig_base()
634 		if (up->maxsignal < DRPOUT)  in irig_base()
635 			up->errflg |= IRIG_ERR_AMP;  in irig_base()
636 		if (up->maxsignal > 0)  in irig_base()
637 			up->modndx = (up->maxsignal - up->noise) /  in irig_base()
638 			    up->maxsignal;  in irig_base()
640 			up->modndx = 0;  in irig_base()
641 		if (up->modndx < MODMIN)  in irig_base()
642 			up->errflg |= IRIG_ERR_MOD;  in irig_base()
643 		if (up->errflg & (IRIG_ERR_AMP | IRIG_ERR_FREQ |  in irig_base()
645 			up->tc = MINTC;  in irig_base()
646 			up->tcount = 0;  in irig_base()
656 		dtemp = up->zxing * up->decim / BAUD;  in irig_base()
657 		up->yxing = dtemp;  in irig_base()
658 		up->zxing = 0.;  in irig_base()
659 		up->phase += dtemp / up->tc;  in irig_base()
660 		up->freq += dtemp / (4. * up->tc * up->tc);  in irig_base()
661 		if (up->freq > MAXFREQ) {  in irig_base()
662 			up->freq = MAXFREQ;  in irig_base()
663 			up->errflg |= IRIG_ERR_FREQ;  in irig_base()
664 		} else if (up->freq < -MAXFREQ) {  in irig_base()
665 			up->freq = -MAXFREQ;  in irig_base()
666 			up->errflg |= IRIG_ERR_FREQ;  in irig_base()
673 	 * negative-going zero crossing at sample 4, the maximum  in irig_base()
678 	 * the slice level and left-shifted in the decoding register.  in irig_base()
683 	lope = (up->lastint[2] - up->lastint[6]) / 2.;  in irig_base()
684 	if (lope > up->intmax)  in irig_base()
685 		up->intmax = lope;  in irig_base()
686 	if (lope < up->intmin)  in irig_base()
687 		up->intmin = lope;  in irig_base()
695 	up->pulse = (up->pulse + 1) % 10;  in irig_base()
696 	up->cycles <<= 1;  in irig_base()
697 	if (lope >= (up->maxsignal + up->noise) / 2.)  in irig_base()
698 		up->cycles |= 1;  in irig_base()
699 	if ((up->cycles & 0x303c0f03) == 0x300c0300) {  in irig_base()
700 		if (up->pulse != 0)  in irig_base()
701 			up->errflg |= IRIG_ERR_SYNCH;  in irig_base()
702 		up->pulse = 0;  in irig_base()
706 	 * Assemble the baud and max/min to get the slice level for the  in irig_base()
707 	 * next baud. The slice level is based on the maximum over the  in irig_base()
709 	 * the slice level halfway between the maximum and minimum.  in irig_base()
711 	env = (up->lastenv[2] - up->lastenv[6]) / 2.;  in irig_base()
712 	up->dcycles <<= 1;  in irig_base()
713 	if (env >= up->slice)  in irig_base()
714 		up->dcycles |= 1;  in irig_base()
715 	switch(up->pulse) {  in irig_base()
718 		irig_baud(peer, up->dcycles);  in irig_base()
719 		if (env < up->envmin)  in irig_base()
720 			up->envmin = env;  in irig_base()
721 		up->slice = (up->envmax + up->envmin) / 2;  in irig_base()
722 		up->envmin = 1e6; up->envmax = -1e6;  in irig_base()
726 		up->envmax = env;  in irig_base()
730 		if (env > up->envmax)  in irig_base()
731 			up->envmax = env;  in irig_base()
735 		up->envmin = env;  in irig_base()
741  * irig_baud - update the PLL and decode the pulse-width signal
754         pp = peer->procptr;  in irig_baud()
755 	up = pp->unitptr;  in irig_baud()
764 	up->exing = -up->yxing;  in irig_baud()
765 	if (abs(up->envxing - up->envphase) <= 1) {  in irig_baud()
766 		up->tcount++;  in irig_baud()
767 		if (up->tcount > 20 * up->tc) {  in irig_baud()
768 			up->tc++;  in irig_baud()
769 			if (up->tc > MAXTC)  in irig_baud()
770 				up->tc = MAXTC;  in irig_baud()
771 			up->tcount = 0;  in irig_baud()
772 			up->envxing = up->envphase;  in irig_baud()
774 			up->exing -= up->envxing - up->envphase;  in irig_baud()
777 		up->tcount = 0;  in irig_baud()
778 		up->envxing = up->envphase;  in irig_baud()
786 	up->prvstamp = up->chrstamp;  in irig_baud()
787 	dtemp = up->decim * (up->exing / SECOND) + up->fdelay;  in irig_baud()
789 	up->chrstamp = up->timestamp;  in irig_baud()
790 	L_SUB(&up->chrstamp, &ltemp);  in irig_baud()
793 	 * The data bits are collected in ten-bit bauds. The first two  in irig_baud()
795 	 * 0-1 bits (0), 2-4 bits (1) and 5-7 bits (PI). The remaining  in irig_baud()
796 	 * 8-bit run represents a soft error and is treated as 0.  in irig_baud()
798 	switch (up->dcycles & 0xff) {  in irig_baud()
800 	case 0x00:		/* 0-1 bits (0) */  in irig_baud()
805 	case 0xc0:		/* 2-4 bits (1) */  in irig_baud()
811 	case 0xf8:		/* (5-7 bits (PI) */  in irig_baud()
819 		up->errflg |= IRIG_ERR_DECODE;  in irig_baud()
825  * irig_decode - decode the data
829  * or position identifier. There are four bits per digit, ten digits per
830  * frame and ten frames per second.
850 	pp = peer->procptr;  in irig_decode()
851 	up = pp->unitptr;  in irig_decode()
856 	up->bits >>= 1;  in irig_decode()
858 		up->bits |= 0x200;  in irig_decode()
859 	} else if (bit == BITP && up->lastbit == BITP) {  in irig_decode()
862 		 * Frame sync - two adjacent position identifiers, which  in irig_decode()
868 		if (up->frmcnt != 1)  in irig_decode()
869 			up->errflg |= IRIG_ERR_SYNCH;  in irig_decode()
870 		up->frmcnt = 1;  in irig_decode()
871 		up->refstamp = up->prvstamp;  in irig_decode()
873 	up->lastbit = bit;  in irig_decode()
874 	if (up->frmcnt % SUBFLD == 0) {  in irig_decode()
878 		 * little-endian timecode field. Note frame 1 is shifted  in irig_decode()
881 		temp = up->bits;  in irig_decode()
882 		if (up->frmcnt == 10)  in irig_decode()
884 		if (up->xptr >= 2) {  in irig_decode()
885 			up->timecode[--up->xptr] = hexchar[temp & 0xf];  in irig_decode()
886 			up->timecode[--up->xptr] = hexchar[(temp >> 5) &  in irig_decode()
889 		if (up->frmcnt == 0) {  in irig_decode()
904 			up->xptr = 2 * SUBFLD;  in irig_decode()
905 			if (sscanf((char *)up->timecode,  in irig_decode()
906 			   "%6s%2d%1d%2s%3d%2d%2d%2d", sbs, &pp->year,  in irig_decode()
907 			    &syncdig, spare, &pp->day, &pp->hour,  in irig_decode()
908 			    &pp->minute, &pp->second) != 8)  in irig_decode()
909 				pp->leap = LEAP_NOTINSYNC;  in irig_decode()
911 				pp->leap = LEAP_NOWARNING;  in irig_decode()
912 			up->second = (up->second + up->decim) % 60;  in irig_decode()
926 			if (pp->day == 0 || (pp->year != 0 && syncdig ==  in irig_decode()
928 				up->errflg |= IRIG_ERR_SIGERR;  in irig_decode()
929 			if (pp->second != up->second)  in irig_decode()
930 				up->errflg |= IRIG_ERR_CHECK;  in irig_decode()
931 			up->second = pp->second;  in irig_decode()
938 			if (up->errflg == 0 && up->tc == MAXTC) {  in irig_decode()
939 				pp->lastref = pp->lastrec;  in irig_decode()
940 				pp->lastrec = up->refstamp;  in irig_decode()
945 			snprintf(pp->a_lastcode, sizeof(pp->a_lastcode),  in irig_decode()
947 			    up->errflg, pp->year, pp->day,  in irig_decode()
948 			    pp->hour, pp->minute, pp->second,  in irig_decode()
949 			    up->maxsignal, up->gain, up->modndx,  in irig_decode()
950 			    up->tc, up->exing * 1e6 / SECOND, up->freq *  in irig_decode()
951 			    1e6 / SECOND, ulfptoa(&pp->lastrec, 6));  in irig_decode()
952 			pp->lencode = strlen(pp->a_lastcode);  in irig_decode()
953 			up->errflg = 0;  in irig_decode()
954 			if (pp->sloppyclockflag & CLK_FLAG4) {  in irig_decode()
955 				record_clock_stats(&peer->srcadr,  in irig_decode()
956 				    pp->a_lastcode);  in irig_decode()
960 					    pp->a_lastcode);  in irig_decode()
965 	up->frmcnt = (up->frmcnt + 1) % FIELD;  in irig_decode()
970  * irig_poll - called by the transmit procedure
973  * IRIG-B there should be at least 60 updates; for IRIG-E there should
984 	pp = peer->procptr;  in irig_poll()
986 	if (pp->coderecv == pp->codeproc) {  in irig_poll()
992 	if (!(pp->sloppyclockflag & CLK_FLAG4)) {  in irig_poll()
993 		record_clock_stats(&peer->srcadr, pp->a_lastcode);  in irig_poll()
996 			printf("irig %s\n", pp->a_lastcode);  in irig_poll()
999 	pp->polls++;  in irig_poll()
1005  * irig_gain - adjust codec gain
1020 	pp = peer->procptr;  in irig_gain()
1021 	up = pp->unitptr;  in irig_gain()
1028 	if (up->maxsignal < MINAMP) {  in irig_gain()
1029 		up->gain += 4;  in irig_gain()
1030 		if (up->gain > MAXGAIN)  in irig_gain()
1031 			up->gain = MAXGAIN;  in irig_gain()
1032 	} else if (up->maxsignal > MAXAMP) {  in irig_gain()
1033 		up->gain -= 4;  in irig_gain()
1034 		if (up->gain < 0)  in irig_gain()
1035 			up->gain = 0;  in irig_gain()
1037 	audio_gain(up->gain, up->mongain, up->port);  in irig_gain()