refclock_irig.c revision 285612
1139823Simp/*
210937Swollman * refclock_irig - audio IRIG-B/E demodulator/decoder
31541Srgrimes */
41541Srgrimes#ifdef HAVE_CONFIG_H
51541Srgrimes#include <config.h>
61541Srgrimes#endif
71541Srgrimes
81541Srgrimes#if defined(REFCLOCK) && defined(CLOCK_IRIG)
91541Srgrimes
101541Srgrimes#include "ntpd.h"
111541Srgrimes#include "ntp_io.h"
121541Srgrimes#include "ntp_refclock.h"
131541Srgrimes#include "ntp_calendar.h"
141541Srgrimes#include "ntp_stdlib.h"
151541Srgrimes
161541Srgrimes#include <stdio.h>
171541Srgrimes#include <ctype.h>
181541Srgrimes#include <math.h>
191541Srgrimes#ifdef HAVE_SYS_IOCTL_H
201541Srgrimes#include <sys/ioctl.h>
211541Srgrimes#endif /* HAVE_SYS_IOCTL_H */
221541Srgrimes
231541Srgrimes#include "audio.h"
241541Srgrimes
251541Srgrimes/*
261541Srgrimes * Audio IRIG-B/E demodulator/decoder
271541Srgrimes *
281541Srgrimes * This driver synchronizes the computer time using data encoded in
2910937Swollman * IRIG-B/E signals commonly produced by GPS receivers and other timing
3050477Speter * devices. The IRIG signal is an amplitude-modulated carrier with
311541Srgrimes * pulse-width modulated data bits. For IRIG-B, the carrier frequency is
321541Srgrimes * 1000 Hz and bit rate 100 b/s; for IRIG-E, the carrier frequenchy is
332169Spaul * 100 Hz and bit rate 10 b/s. The driver automatically recognizes which
342169Spaul & format is in use.
3586764Sjlemon *
3698102Shsu * The driver requires an audio codec or sound card with sampling rate 8
3786764Sjlemon * kHz and mu-law companding. This is the same standard as used by the
381541Srgrimes * telephone industry and is supported by most hardware and operating
391541Srgrimes * systems, including Solaris, SunOS, FreeBSD, NetBSD and Linux. In this
401541Srgrimes * implementation, only one audio driver and codec can be supported on a
4186764Sjlemon * single machine.
421541Srgrimes *
4355679Sshin * The program processes 8000-Hz mu-law companded samples using separate
4455679Sshin * signal filters for IRIG-B and IRIG-E, a comb filter, envelope
4560938Sjake * detector and automatic threshold corrector. Cycle crossings relative
4655679Sshin * to the corrected slice level determine the width of each pulse and
4755679Sshin * its value - zero, one or position identifier.
4855679Sshin *
4955679Sshin * The data encode 20 BCD digits which determine the second, minute,
5060938Sjake * hour and day of the year and sometimes the year and synchronization
51126193Sandre * condition. The comb filter exponentially averages the corresponding
52126225Sbde * samples of successive baud intervals in order to reliably identify
5355679Sshin * the reference carrier cycle. A type-II phase-lock loop (PLL) performs
54130989Sps * additional integration and interpolation to accurately determine the
55130989Sps * zero crossing of that cycle, which determines the reference
56131078Sbms * timestamp. A pulse-width discriminator demodulates the data pulses,
57130989Sps * which are then encoded as the BCD digits of the timecode.
58130989Sps *
59130989Sps * The timecode and reference timestamp are updated once each second
60130989Sps * with IRIG-B (ten seconds with IRIG-E) and local clock offset samples
61130989Sps * saved for later processing. At poll intervals of 64 s, the saved
62130989Sps * samples are processed by a trimmed-mean filter and used to update the
63145370Sps * system clock.
64130989Sps *
65131078Sbms * An automatic gain control feature provides protection against
66146123Sps * overdriven or underdriven input signal amplitudes. It is designed to
67146123Sps * maintain adequate demodulator signal amplitude while avoiding
68146123Sps * occasional noise spikes. In order to assure reliable capture, the
69146123Sps * decompanded input signal amplitude must be greater than 100 units and
70146123Sps * the codec sample frequency error less than 250 PPM (.025 percent).
7155679Sshin *
7255679Sshin * Monitor Data
7355679Sshin *
7455679Sshin * The timecode format used for debugging and data recording includes
7555679Sshin * data helpful in diagnosing problems with the IRIG signal and codec
7655679Sshin * connections. The driver produces one line for each timecode in the
7755679Sshin * following format:
78169541Sandre *
79169541Sandre * 00 00 98 23 19:26:52 2782 143 0.694 10 0.3 66.5 3094572411.00027
80169541Sandre *
81169541Sandre * If clockstats is enabled, the most recent line is written to the
82169541Sandre * clockstats file every 64 s. If verbose recording is enabled (fudge
83169541Sandre * flag 4) each line is written as generated.
84169541Sandre *
85169541Sandre * The first field containes the error flags in hex, where the hex bits
86169541Sandre * are interpreted as below. This is followed by the year of century,
87169541Sandre * day of year and time of day. Note that the time of day is for the
88169541Sandre * previous minute, not the current time. The status indicator and year
89169541Sandre * are not produced by some IRIG devices and appear as zeros. Following
901541Srgrimes * these fields are the carrier amplitude (0-3000), codec gain (0-255),
911541Srgrimes * modulation index (0-1), time constant (4-10), carrier phase error
9232821Sdg * +-.5) and carrier frequency error (PPM). The last field is the on-
931541Srgrimes * time timestamp in NTP format.
941541Srgrimes *
95126193Sandre * The error flags are defined as follows in hex:
96126193Sandre *
9732821Sdg * x01	Low signal. The carrier amplitude is less than 100 units. This
9832821Sdg *	is usually the result of no signal or wrong input port.
99168615Sandre * x02	Frequency error. The codec frequency error is greater than 250
10032821Sdg *	PPM. This may be due to wrong signal format or (rarely)
10132821Sdg *	defective codec.
10211187Swollman * x04	Modulation error. The IRIG modulation index is less than 0.5.
10311187Swollman *	This is usually the result of an overdriven codec, wrong signal
104117650Shsu *	format or wrong input port.
105117650Shsu * x08	Frame synch error. The decoder frame does not match the IRIG
106117650Shsu *	frame. This is usually the result of an overdriven codec, wrong
107117650Shsu *	signal format or noisy IRIG signal. It may also be the result of
108117650Shsu *	an IRIG signature check which indicates a failure of the IRIG
109117650Shsu *	signal synchronization source.
110117650Shsu * x10	Data bit error. The data bit length is out of tolerance. This is
111117650Shsu *	usually the result of an overdriven codec, wrong signal format
112117650Shsu *	or noisy IRIG signal.
113117650Shsu * x20	Seconds numbering discrepancy. The decoder second does not match
114117650Shsu *	the IRIG second. This is usually the result of an overdriven
115117650Shsu *	codec, wrong signal format or noisy IRIG signal.
116117650Shsu * x40	Codec error (overrun). The machine is not fast enough to keep up
117117650Shsu *	with the codec.
118117650Shsu * x80	Device status error (Spectracom).
119117650Shsu *
120117650Shsu *
121117650Shsu * Once upon a time, an UltrSPARC 30 and Solaris 2.7 kept the clock
122117650Shsu * within a few tens of microseconds relative to the IRIG-B signal.
123125680Sbms * Accuracy with IRIG-E was about ten times worse. Unfortunately, Sun
124146463Sps * broke the 2.7 audio driver in 2.8, which has a 10-ms sawtooth
125162084Sandre * modulation.
1261541Srgrimes *
1271541Srgrimes * Unlike other drivers, which can have multiple instantiations, this
12832821Sdg * one supports only one. It does not seem likely that more than one
12932821Sdg * audio codec would be useful in a single machine. More than one would
13032821Sdg * probably chew up too much CPU time anyway.
1311541Srgrimes *
1321541Srgrimes * Fudge factors
13332821Sdg *
1341541Srgrimes * Fudge flag4 causes the dubugging output described above to be
1351541Srgrimes * recorded in the clockstats file. Fudge flag2 selects the audio input
1361541Srgrimes * port, where 0 is the mike port (default) and 1 is the line-in port.
13732821Sdg * It does not seem useful to select the compact disc player port. Fudge
13832821Sdg * flag3 enables audio monitoring of the input signal. For this purpose,
13932821Sdg * the monitor gain is set t a default value. Fudgetime2 is used as a
14032821Sdg * frequency vernier for broken codec sample frequency.
1411541Srgrimes *
1421541Srgrimes * Alarm codes
14332821Sdg *
14432821Sdg * CEVNT_BADTIME	invalid date or time
1451541Srgrimes * CEVNT_TIMEOUT	no IRIG data since last poll
146102017Sdillon */
14713765Smpp/*
1481541Srgrimes * Interface definitions
1491541Srgrimes */
1501541Srgrimes#define	DEVICE_AUDIO	"/dev/audio" /* audio device name */
151102017Sdillon#define	PRECISION	(-17)	/* precision assumed (about 10 us) */
152109175Shsu#define	REFID		"IRIG"	/* reference ID */
15360067Sjlemon#define	DESCRIPTION	"Generic IRIG Audio Driver" /* WRU */
15432821Sdg#define	AUDIO_BUFSIZ	320	/* audio buffer size (40 ms) */
15532821Sdg#define SECOND		8000	/* nominal sample rate (Hz) */
15650673Sjlemon#define BAUD		80	/* samples per baud interval */
15750673Sjlemon#define OFFSET		128	/* companded sample offset */
15850673Sjlemon#define SIZE		256	/* decompanding table size */
1591541Srgrimes#define CYCLE		8	/* samples per bit */
16032821Sdg#define SUBFLD		10	/* bits per frame */
161102017Sdillon#define FIELD		100	/* bits per second */
162102017Sdillon#define MINTC		2	/* min PLL time constant */
163102017Sdillon#define MAXTC		10	/* max PLL time constant max */
16450673Sjlemon#define	MAXAMP		3000.	/* maximum signal amplitude */
16532821Sdg#define	MINAMP		2000.	/* minimum signal amplitude */
16611187Swollman#define DRPOUT		100.	/* dropout signal amplitude */
16711187Swollman#define MODMIN		0.5	/* minimum modulation index */
16832821Sdg#define MAXFREQ		(250e-6 * SECOND) /* freq tolerance (.025%) */
16932821Sdg
17011187Swollman/*
171102017Sdillon * The on-time synchronization point is the positive-going zero crossing
17232821Sdg * of the first cycle of the second. The IIR baseband filter phase delay
1731541Srgrimes * is 1.03 ms for IRIG-B and 3.47 ms for IRIG-E. The fudge value 2.68 ms
1741541Srgrimes * due to the codec and other causes was determined by calibrating to a
17532821Sdg * PPS signal from a GPS receiver.
1761541Srgrimes *
1771541Srgrimes * The results with a 2.4-GHz P4 running FreeBSD 6.1 are generally
1781541Srgrimes * within .02 ms short-term with .02 ms jitter. The processor load due
1791541Srgrimes * to the driver is 0.51 percent.
1801541Srgrimes */
1811541Srgrimes#define IRIG_B	((1.03 + 2.68) / 1000)	/* IRIG-B system delay (s) */
1821541Srgrimes#define IRIG_E	((3.47 + 2.68) / 1000)	/* IRIG-E system delay (s) */
1831541Srgrimes
1841541Srgrimes/*
185162277Sandre * Data bit definitions
186162277Sandre */
187162277Sandre#define BIT0		0	/* zero */
18832821Sdg#define BIT1		1	/* one */
1891541Srgrimes#define BITP		2	/* position identifier */
19050673Sjlemon
19150673Sjlemon/*
19250673Sjlemon * Error flags
193117650Shsu */
19450673Sjlemon#define IRIG_ERR_AMP	0x01	/* low carrier amplitude */
195112957Shsu#define IRIG_ERR_FREQ	0x02	/* frequency tolerance exceeded */
196131078Sbms#define IRIG_ERR_MOD	0x04	/* low modulation index */
197130989Sps#define IRIG_ERR_SYNCH	0x08	/* frame synch error */
198146953Sps#define IRIG_ERR_DECODE	0x10	/* frame decoding error */
199146953Sps#define IRIG_ERR_CHECK	0x20	/* second numbering discrepancy */
200146630Sps#define IRIG_ERR_ERROR	0x40	/* codec error (overrun) */
201130989Sps#define IRIG_ERR_SIGERR	0x80	/* IRIG status error (Spectracom) */
202130989Sps
203136151Spsstatic	char	hexchar[] = "0123456789abcdef";
204136151Sps
205146123Sps/*
206155767Sandre * IRIG unit control structure
207166405Sandre */
208166405Sandrestruct irigunit {
2091541Srgrimes	u_char	timecode[2 * SUBFLD + 1]; /* timecode string */
2101541Srgrimes	l_fp	timestamp;	/* audio sample timestamp */
211117650Shsu	l_fp	tick;		/* audio sample increment */
212117650Shsu	l_fp	refstamp;	/* reference timestamp */
213117650Shsu	l_fp	chrstamp;	/* baud timestamp */
214117650Shsu	l_fp	prvstamp;	/* previous baud timestamp */
215125680Sbms	double	integ[BAUD];	/* baud integrator */
2166247Swollman	double	phase, freq;	/* logical clock phase and frequency */
217125680Sbms	double	zxing;		/* phase detector integrator */
218125680Sbms	double	yxing;		/* cycle phase */
219125680Sbms	double	exing;		/* envelope phase */
220125680Sbms	double	modndx;		/* modulation index */
221125680Sbms	double	irig_b;		/* IRIG-B signal amplitude */
222125680Sbms	double	irig_e;		/* IRIG-E signal amplitude */
223125680Sbms	int	errflg;		/* error flags */
224125680Sbms	/*
225125680Sbms	 * Audio codec variables
226125680Sbms	 */
227125680Sbms	double	comp[SIZE];	/* decompanding table */
228125680Sbms	double	signal;		/* peak signal for AGC */
229125680Sbms	int	port;		/* codec port */
230125680Sbms	int	gain;		/* codec gain */
2316247Swollman	int	mongain;	/* codec monitor gain */
2326247Swollman	int	seccnt;		/* second interval counter */
2336247Swollman
234167606Sandre	/*
235167606Sandre	 * RF variables
236167606Sandre	 */
2376247Swollman	double	bpf[9];		/* IRIG-B filter shift register */
2386247Swollman	double	lpf[5];		/* IRIG-E filter shift register */
23986764Sjlemon	double	envmin, envmax;	/* envelope min and max */
240167606Sandre	double	slice;		/* envelope slice level */
241167606Sandre	double	intmin, intmax;	/* integrated envelope min and max */
242167606Sandre	double	maxsignal;	/* integrated peak amplitude */
243167606Sandre	double	noise;		/* integrated noise amplitude */
244168906Sandre	double	lastenv[CYCLE];	/* last cycle amplitudes */
245168906Sandre	double	lastint[CYCLE];	/* last integrated cycle amplitudes */
246168906Sandre	double	lastsig;	/* last carrier sample */
247168906Sandre	double	fdelay;		/* filter delay */
248167606Sandre	int	decim;		/* sample decimation factor */
249167606Sandre	int	envphase;	/* envelope phase */
250167606Sandre	int	envptr;		/* envelope phase pointer */
251147637Sps	int	envsw;		/* envelope state */
252147637Sps	int	envxing;	/* envelope slice crossing */
253168906Sandre	int	tc;		/* time constant */
2546247Swollman	int	tcount;		/* time constant counter */
2556247Swollman	int	badcnt;		/* decimation interval counter */
256159949Sandre
257159949Sandre	/*
258159949Sandre	 * Decoder variables
259159949Sandre	 */
260159949Sandre	int	pulse;		/* cycle counter */
261122922Sandre	int	cycles;		/* carrier cycles */
262122922Sandre	int	dcycles;	/* data cycles */
263122922Sandre	int	lastbit;	/* last code element */
264122922Sandre	int	second;		/* previous second */
265122922Sandre	int	bitcnt;		/* bit count in frame */
266122922Sandre	int	frmcnt;		/* bit count in second */
267122922Sandre	int	xptr;		/* timecode pointer */
268122922Sandre	int	bits;		/* demodulated bits */
269122922Sandre};
270122922Sandre
271122922Sandre/*
272159725Sandre * Function prototypes
273159725Sandre */
274159725Sandrestatic	int	irig_start	(int, struct peer *);
275159725Sandrestatic	void	irig_shutdown	(int, struct peer *);
276111145Sjlemonstatic	void	irig_receive	(struct recvbuf *);
277111145Sjlemonstatic	void	irig_poll	(int, struct peer *);
278111145Sjlemon
279111145Sjlemon/*
280121850Ssilby * More function prototypes
281121884Ssilby */
282111145Sjlemonstatic	void	irig_base	(struct peer *, double);
283111145Sjlemonstatic	void	irig_rf		(struct peer *, double);
284111145Sjlemonstatic	void	irig_baud	(struct peer *, int);
285133874Srwatsonstatic	void	irig_decode	(struct peer *, int);
286169477Sandrestatic	void	irig_gain	(struct peer *);
287111145Sjlemon
288112009Sjlemon/*
289162111Sru * Transfer vector
290111145Sjlemon */
291133874Srwatsonstruct	refclock refclock_irig = {
2921541Srgrimes	irig_start,		/* start up driver */
293111145Sjlemon	irig_shutdown,		/* shut down driver */
2941541Srgrimes	irig_poll,		/* transmit poll message */
2951541Srgrimes	noentry,		/* not used (old irig_control) */
2961541Srgrimes	noentry,		/* initialize driver (not used) */
2971541Srgrimes	noentry,		/* not used (old irig_buginfo) */
2981541Srgrimes	NOFLAGS			/* not used */
2991541Srgrimes};
3001541Srgrimes
3011541Srgrimes
3021541Srgrimes/*
3031541Srgrimes * irig_start - open the devices and initialize data for processing
3041541Srgrimes */
30514753Swollmanstatic int
30614753Swollmanirig_start(
30714753Swollman	int	unit,		/* instance number (used for PCM) */
30814753Swollman	struct peer *peer	/* peer structure pointer */
30914753Swollman	)
3101541Srgrimes{
3111541Srgrimes	struct refclockproc *pp;
3121541Srgrimes	struct irigunit *up;
3131541Srgrimes
3141541Srgrimes	/*
3151541Srgrimes	 * Local variables
3161541Srgrimes	 */
3171541Srgrimes	int	fd;		/* file descriptor */
3181541Srgrimes	int	i;		/* index */
3191541Srgrimes	double	step;		/* codec adjustment */
3201541Srgrimes
32114753Swollman	/*
32214753Swollman	 * Open audio device
32314753Swollman	 */
32414753Swollman	fd = audio_init(DEVICE_AUDIO, AUDIO_BUFSIZ, unit);
32514753Swollman	if (fd < 0)
3261541Srgrimes		return (0);
3271541Srgrimes#ifdef DEBUG
32835419Sdg	if (debug)
32916141Swollman		audio_show();
3301541Srgrimes#endif
3311541Srgrimes
3321541Srgrimes	/*
3331541Srgrimes	 * Allocate and initialize unit structure
3341541Srgrimes	 */
3351541Srgrimes	up = emalloc_zero(sizeof(*up));
3361541Srgrimes	pp = peer->procptr;
3371541Srgrimes	pp->io.clock_recv = irig_receive;
3381541Srgrimes	pp->io.srcclock = peer;
3391541Srgrimes	pp->io.datalen = 0;
3401541Srgrimes	pp->io.fd = fd;
3411541Srgrimes	if (!io_addclock(&pp->io)) {
342124258Sandre		close(fd);
3431541Srgrimes		pp->io.fd = -1;
3441541Srgrimes		free(up);
3451541Srgrimes		return (0);
3461541Srgrimes	}
3471541Srgrimes	pp->unitptr = up;
3481541Srgrimes
3491541Srgrimes	/*
3501541Srgrimes	 * Initialize miscellaneous variables
3511541Srgrimes	 */
3521541Srgrimes	peer->precision = PRECISION;
3531541Srgrimes	pp->clockdesc = DESCRIPTION;
3541541Srgrimes	memcpy((char *)&pp->refid, REFID, 4);
3551541Srgrimes	up->tc = MINTC;
3561541Srgrimes	up->decim = 1;
3571541Srgrimes	up->gain = 127;
3581541Srgrimes
359100373Sdillon	/*
3601541Srgrimes	 * The companded samples are encoded sign-magnitude. The table
3611541Srgrimes	 * contains all the 256 values in the interest of speed.
3621541Srgrimes	 */
3631541Srgrimes	up->comp[0] = up->comp[OFFSET] = 0.;
3641541Srgrimes	up->comp[1] = 1; up->comp[OFFSET + 1] = -1.;
3651541Srgrimes	up->comp[2] = 3; up->comp[OFFSET + 2] = -3.;
3661541Srgrimes	step = 2.;
3671541Srgrimes	for (i = 3; i < OFFSET; i++) {
3681541Srgrimes		up->comp[i] = up->comp[i - 1] + step;
3691541Srgrimes		up->comp[OFFSET + i] = -up->comp[i];
3701541Srgrimes		if (i % 16 == 0)
37152904Sshin			step *= 2.;
3721541Srgrimes	}
3731541Srgrimes	DTOLFP(1. / SECOND, &up->tick);
3741541Srgrimes	return (1);
3751541Srgrimes}
3761541Srgrimes
3771541Srgrimes
3781541Srgrimes/*
3791541Srgrimes * irig_shutdown - shut down the clock
3801541Srgrimes */
3811541Srgrimesstatic void
3821541Srgrimesirig_shutdown(
3831541Srgrimes	int	unit,		/* instance number (not used) */
3841541Srgrimes	struct peer *peer	/* peer structure pointer */
3851541Srgrimes	)
3861541Srgrimes{
3871541Srgrimes	struct refclockproc *pp;
3881541Srgrimes	struct irigunit *up;
3891541Srgrimes
3901541Srgrimes	pp = peer->procptr;
3911541Srgrimes	up = pp->unitptr;
3929263Swollman	if (-1 != pp->io.fd)
3939263Swollman		io_closeclock(&pp->io);
3949263Swollman	if (NULL != up)
3959470Swollman		free(up);
3969470Swollman}
3979470Swollman
39810937Swollman
39910937Swollman/*
40011415Swollman * irig_receive - receive data from the audio device
40114281Sbde *
402128653Ssilby * This routine reads input samples and adjusts the logical clock to
40386764Sjlemon * track the irig clock by dropping or duplicating codec samples.
40486764Sjlemon */
40586764Sjlemonstatic void
40686764Sjlemonirig_receive(
40786764Sjlemon	struct recvbuf *rbufp	/* receive buffer structure pointer */
40886764Sjlemon	)
40986764Sjlemon{
41086764Sjlemon	struct peer *peer;
41186764Sjlemon	struct refclockproc *pp;
41286764Sjlemon	struct irigunit *up;
41386764Sjlemon
41486764Sjlemon	/*
41586764Sjlemon	 * Local variables
41686764Sjlemon	 */
41786764Sjlemon	double	sample;		/* codec sample */
41886764Sjlemon	u_char	*dpt;		/* buffer pointer */
419122922Sandre	int	bufcnt;		/* buffer counter */
420122922Sandre	l_fp	ltemp;		/* l_fp temp */
421122922Sandre
422130989Sps	peer = rbufp->recv_peer;
423167036Smohans	pp = peer->procptr;
424167036Smohans	up = pp->unitptr;
425130989Sps
426130989Sps	/*
427130989Sps	 * Main loop - read until there ain't no more. Note codec
428133874Srwatson	 * samples are bit-inverted.
429130989Sps	 */
430130989Sps	DTOLFP((double)rbufp->recv_length / SECOND, &ltemp);
431143339Sps	L_SUB(&rbufp->recv_time, &ltemp);
4321541Srgrimes	up->timestamp = rbufp->recv_time;
4331541Srgrimes	dpt = rbufp->recv_buffer;
4346247Swollman	for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) {
43536079Swollman		sample = up->comp[~*dpt++ & 0xff];
43637183Sjhay
43737183Sjhay		/*
43836079Swollman		 * Variable frequency oscillator. The codec oscillator
43937183Sjhay		 * runs at the nominal rate of 8000 samples per second,
44036079Swollman		 * or 125 us per sample. A frequency change of one unit
44136079Swollman		 * results in either duplicating or deleting one sample
44236079Swollman		 * per second, which results in a frequency change of
44336079Swollman		 * 125 PPM.
44436079Swollman		 */
44536079Swollman		up->phase += (up->freq + clock_codec) / SECOND;
44636079Swollman		up->phase += pp->fudgetime2 / 1e6;
44736079Swollman		if (up->phase >= .5) {
44836079Swollman			up->phase -= 1.;
44936079Swollman		} else if (up->phase < -.5) {
4506247Swollman			up->phase += 1.;
4516247Swollman			irig_rf(peer, sample);
4526247Swollman			irig_rf(peer, sample);
4536247Swollman		} else {
4546472Swollman			irig_rf(peer, sample);
4556472Swollman		}
4566472Swollman		L_ADD(&up->timestamp, &up->tick);
4576472Swollman		sample = fabs(sample);
45818281Spst		if (sample > up->signal)
45918281Spst			up->signal = sample;
46063431Ssheldonh		up->signal += (sample - up->signal) /
46136079Swollman		    1000;
46250673Sjlemon
46352904Sshin		/*
464130989Sps		 * Once each second, determine the IRIG format and gain.
465141381Smaxim		 */
466141381Smaxim		up->seccnt = (up->seccnt + 1) % SECOND;
467167036Smohans		if (up->seccnt == 0) {
4686247Swollman			if (up->irig_b > up->irig_e) {
4696247Swollman				up->decim = 1;
4706247Swollman				up->fdelay = IRIG_B;
4716348Swollman			} else {
4726247Swollman				up->decim = 10;
4736472Swollman				up->fdelay = IRIG_E;
4746472Swollman			}
4756472Swollman			up->irig_b = up->irig_e = 0;
4766472Swollman			irig_gain(peer);
4776472Swollman
4786472Swollman		}
47918281Spst	}
48036079Swollman
48150673Sjlemon	/*
48252904Sshin	 * Set the input port and monitor gain for the next buffer.
483122922Sandre	 */
4846247Swollman	if (pp->sloppyclockflag & CLK_FLAG2)
4856247Swollman		up->port = 2;
48644078Sdfr	else
48755205Speter		up->port = 1;
48844078Sdfr	if (pp->sloppyclockflag & CLK_FLAG3)
48944078Sdfr		up->mongain = MONGAIN;
490136151Sps	else
49144078Sdfr		up->mongain = 0;
49244078Sdfr}
4937684Sdg
4947684Sdg
4956348Swollman/*
4967090Sbde * irig_rf - RF processing
497124258Sandre *
49833846Sdg * This routine filters the RF signal using a bandass filter for IRIG-B
49960067Sjlemon * and a lowpass filter for IRIG-E. In case of IRIG-E, the samples are
50098204Ssilby * decimated by a factor of ten. Note that the codec filters function as
50150673Sjlemon * roofing filters to attenuate both the high and low ends of the
50250673Sjlemon * passband. IIR filter coefficients were determined using Matlab Signal
5031541Srgrimes * Processing Toolkit.
504130989Sps */
505130989Spsstatic void
506167606Sandreirig_rf(
5071541Srgrimes	struct peer *peer,	/* peer structure pointer */
50892723Salfred	double	sample		/* current signal sample */
509157376Srwatson	)
510111145Sjlemon{
511162064Sglebius	struct refclockproc *pp;
512121850Ssilby	struct irigunit *up;
513162064Sglebius
514157376Srwatson	/*
51592723Salfred	 * Local variables
51692723Salfred	 */
5171541Srgrimes	double	irig_b, irig_e;	/* irig filter outputs */
51892723Salfred
51992723Salfred	pp = peer->procptr;
52092723Salfred	up = pp->unitptr;
52192723Salfred
522128452Ssilby	/*
523169541Sandre	 * IRIG-B filter. Matlab 4th-order IIR elliptic, 800-1200 Hz
524126193Sandre	 * bandpass, 0.3 dB passband ripple, -50 dB stopband ripple,
52592723Salfred	 * phase delay 1.03 ms.
526162084Sandre	 */
527162084Sandre	irig_b = (up->bpf[8] = up->bpf[7]) * 6.505491e-001;
52892723Salfred	irig_b += (up->bpf[7] = up->bpf[6]) * -3.875180e+000;
529122922Sandre	irig_b += (up->bpf[6] = up->bpf[5]) * 1.151180e+001;
530133874Srwatson	irig_b += (up->bpf[5] = up->bpf[4]) * -2.141264e+001;
53198211Shsu	irig_b += (up->bpf[4] = up->bpf[3]) * 2.712837e+001;
53298211Shsu	irig_b += (up->bpf[3] = up->bpf[2]) * -2.384486e+001;
53398211Shsu	irig_b += (up->bpf[2] = up->bpf[1]) * 1.427663e+001;
5341541Srgrimes	irig_b += (up->bpf[1] = up->bpf[0]) * -5.352734e+000;
53592723Salfred	up->bpf[0] = sample - irig_b;
53692723Salfred	irig_b = up->bpf[0] * 4.952157e-003
53792723Salfred	    + up->bpf[1] * -2.055878e-002
53892723Salfred	    + up->bpf[2] * 4.401413e-002
539169541Sandre	    + up->bpf[3] * -6.558851e-002
540169541Sandre	    + up->bpf[4] * 7.462108e-002
541126351Srwatson	    + up->bpf[5] * -6.558851e-002
54292723Salfred	    + up->bpf[6] * 4.401413e-002
543125783Sbms	    + up->bpf[7] * -2.055878e-002
544125783Sbms	    + up->bpf[8] * 4.952157e-003;
545125783Sbms	up->irig_b += irig_b * irig_b;
54692723Salfred
54755679Sshin	/*
548111144Sjlemon	 * IRIG-E filter. Matlab 4th-order IIR elliptic, 130-Hz lowpass,
549111144Sjlemon	 * 0.3 dB passband ripple, -50 dB stopband ripple, phase delay
550168615Sandre	 * 3.47 ms.
551168615Sandre	 */
552169272Srwatson	irig_e = (up->lpf[4] = up->lpf[3]) * 8.694604e-001;
553102017Sdillon	irig_e += (up->lpf[3] = up->lpf[2]) * -3.589893e+000;
55486764Sjlemon	irig_e += (up->lpf[2] = up->lpf[1]) * 5.570154e+000;
55586764Sjlemon	irig_e += (up->lpf[1] = up->lpf[0]) * -3.849667e+000;
556162277Sandre	up->lpf[0] = sample - irig_e;
557162277Sandre	irig_e = up->lpf[0] * 3.215696e-003
558168903Sandre	    + up->lpf[1] * -1.174951e-002
559159695Sandre	    + up->lpf[2] * 1.712074e-002
56086764Sjlemon	    + up->lpf[3] * -1.174951e-002
56186764Sjlemon	    + up->lpf[4] * 3.215696e-003;
562122922Sandre	up->irig_e += irig_e * irig_e;
563122922Sandre
564122922Sandre	/*
565122922Sandre	 * Decimate by a factor of either 1 (IRIG-B) or 10 (IRIG-E).
566122922Sandre	 */
567122922Sandre	up->badcnt = (up->badcnt + 1) % up->decim;
568122922Sandre	if (up->badcnt == 0) {
569122922Sandre		if (up->decim == 1)
5709470Swollman			irig_base(peer, irig_b);
57117096Swollman		else
5729470Swollman			irig_base(peer, irig_e);
5739470Swollman	}
57492723Salfred}
5759470Swollman
576147637Sps/*
577142031Sps * irig_base - baseband processing
578130989Sps *
579130989Sps * This routine processes the baseband signal and demodulates the AM
580136151Sps * carrier using a synchronous detector. It then synchronizes to the
581130989Sps * data frame at the baud rate and decodes the width-modulated data
582130989Sps * pulses.
583130989Sps */
584130989Spsstatic void
585130989Spsirig_base(
58655205Speter	struct peer *peer,	/* peer structure pointer */
5872169Spaul	double	sample		/* current signal sample */
5886348Swollman	)
589{
590	struct refclockproc *pp;
591	struct irigunit *up;
592
593	/*
594	 * Local variables
595	 */
596	double	lope;		/* integrator output */
597	double	env;		/* envelope detector output */
598	double	dtemp;
599	int	carphase;	/* carrier phase */
600
601	pp = peer->procptr;
602	up = pp->unitptr;
603
604	/*
605	 * Synchronous baud integrator. Corresponding samples of current
606	 * and past baud intervals are integrated to refine the envelope
607	 * amplitude and phase estimate. We keep one cycle (1 ms) of the
608	 * raw data and one baud (10 ms) of the integrated data.
609	 */
610	up->envphase = (up->envphase + 1) % BAUD;
611	up->integ[up->envphase] += (sample - up->integ[up->envphase]) /
612	    (5 * up->tc);
613	lope = up->integ[up->envphase];
614	carphase = up->envphase % CYCLE;
615	up->lastenv[carphase] = sample;
616	up->lastint[carphase] = lope;
617
618	/*
619	 * Phase detector. Find the negative-going zero crossing
620	 * relative to sample 4 in the 8-sample sycle. A phase change of
621	 * 360 degrees produces an output change of one unit.
622	 */
623	if (up->lastsig > 0 && lope <= 0)
624		up->zxing += (double)(carphase - 4) / CYCLE;
625	up->lastsig = lope;
626
627	/*
628	 * End of the baud. Update signal/noise estimates and PLL
629	 * phase, frequency and time constant.
630	 */
631	if (up->envphase == 0) {
632		up->maxsignal = up->intmax; up->noise = up->intmin;
633		up->intmin = 1e6; up->intmax = -1e6;
634		if (up->maxsignal < DRPOUT)
635			up->errflg |= IRIG_ERR_AMP;
636		if (up->maxsignal > 0)
637			up->modndx = (up->maxsignal - up->noise) /
638			    up->maxsignal;
639 		else
640			up->modndx = 0;
641		if (up->modndx < MODMIN)
642			up->errflg |= IRIG_ERR_MOD;
643		if (up->errflg & (IRIG_ERR_AMP | IRIG_ERR_FREQ |
644		   IRIG_ERR_MOD | IRIG_ERR_SYNCH)) {
645			up->tc = MINTC;
646			up->tcount = 0;
647		}
648
649		/*
650		 * Update PLL phase and frequency. The PLL time constant
651		 * is set initially to stabilize the frequency within a
652		 * minute or two, then increases to the maximum. The
653		 * frequency is clamped so that the PLL capture range
654		 * cannot be exceeded.
655		 */
656		dtemp = up->zxing * up->decim / BAUD;
657		up->yxing = dtemp;
658		up->zxing = 0.;
659		up->phase += dtemp / up->tc;
660		up->freq += dtemp / (4. * up->tc * up->tc);
661		if (up->freq > MAXFREQ) {
662			up->freq = MAXFREQ;
663			up->errflg |= IRIG_ERR_FREQ;
664		} else if (up->freq < -MAXFREQ) {
665			up->freq = -MAXFREQ;
666			up->errflg |= IRIG_ERR_FREQ;
667		}
668	}
669
670	/*
671	 * Synchronous demodulator. There are eight samples in the cycle
672	 * and ten cycles in the baud. Since the PLL has aligned the
673	 * negative-going zero crossing at sample 4, the maximum
674	 * amplitude is at sample 2 and minimum at sample 6. The
675	 * beginning of the data pulse is determined from the integrated
676	 * samples, while the end of the pulse is determined from the
677	 * raw samples. The raw data bits are demodulated relative to
678	 * the slice level and left-shifted in the decoding register.
679	 */
680	if (carphase != 7)
681		return;
682
683	lope = (up->lastint[2] - up->lastint[6]) / 2.;
684	if (lope > up->intmax)
685		up->intmax = lope;
686	if (lope < up->intmin)
687		up->intmin = lope;
688
689	/*
690	 * Pulse code demodulator and reference timestamp. The decoder
691	 * looks for a sequence of ten bits; the first two bits must be
692	 * one, the last two bits must be zero. Frame synch is asserted
693	 * when three correct frames have been found.
694	 */
695	up->pulse = (up->pulse + 1) % 10;
696	up->cycles <<= 1;
697	if (lope >= (up->maxsignal + up->noise) / 2.)
698		up->cycles |= 1;
699	if ((up->cycles & 0x303c0f03) == 0x300c0300) {
700		if (up->pulse != 0)
701			up->errflg |= IRIG_ERR_SYNCH;
702		up->pulse = 0;
703	}
704
705	/*
706	 * Assemble the baud and max/min to get the slice level for the
707	 * next baud. The slice level is based on the maximum over the
708	 * first two bits and the minimum over the last two bits, with
709	 * the slice level halfway between the maximum and minimum.
710	 */
711	env = (up->lastenv[2] - up->lastenv[6]) / 2.;
712	up->dcycles <<= 1;
713	if (env >= up->slice)
714		up->dcycles |= 1;
715	switch(up->pulse) {
716
717	case 0:
718		irig_baud(peer, up->dcycles);
719		if (env < up->envmin)
720			up->envmin = env;
721		up->slice = (up->envmax + up->envmin) / 2;
722		up->envmin = 1e6; up->envmax = -1e6;
723		break;
724
725	case 1:
726		up->envmax = env;
727		break;
728
729	case 2:
730		if (env > up->envmax)
731			up->envmax = env;
732		break;
733
734	case 9:
735		up->envmin = env;
736		break;
737	}
738}
739
740/*
741 * irig_baud - update the PLL and decode the pulse-width signal
742 */
743static void
744irig_baud(
745	struct peer *peer,	/* peer structure pointer */
746	int	bits		/* decoded bits */
747	)
748{
749	struct refclockproc *pp;
750	struct irigunit *up;
751	double	dtemp;
752	l_fp	ltemp;
753
754        pp = peer->procptr;
755	up = pp->unitptr;
756
757	/*
758	 * The PLL time constant starts out small, in order to
759	 * sustain a frequency tolerance of 250 PPM. It
760	 * gradually increases as the loop settles down. Note
761	 * that small wiggles are not believed, unless they
762	 * persist for lots of samples.
763	 */
764	up->exing = -up->yxing;
765	if (abs(up->envxing - up->envphase) <= 1) {
766		up->tcount++;
767		if (up->tcount > 20 * up->tc) {
768			up->tc++;
769			if (up->tc > MAXTC)
770				up->tc = MAXTC;
771			up->tcount = 0;
772			up->envxing = up->envphase;
773		} else {
774			up->exing -= up->envxing - up->envphase;
775		}
776	} else {
777		up->tcount = 0;
778		up->envxing = up->envphase;
779	}
780
781	/*
782	 * Strike the baud timestamp as the positive zero crossing of
783	 * the first bit, accounting for the codec delay and filter
784	 * delay.
785	 */
786	up->prvstamp = up->chrstamp;
787	dtemp = up->decim * (up->exing / SECOND) + up->fdelay;
788	DTOLFP(dtemp, &ltemp);
789	up->chrstamp = up->timestamp;
790	L_SUB(&up->chrstamp, &ltemp);
791
792	/*
793	 * The data bits are collected in ten-bit bauds. The first two
794	 * bits are not used. The resulting patterns represent runs of
795	 * 0-1 bits (0), 2-4 bits (1) and 5-7 bits (PI). The remaining
796	 * 8-bit run represents a soft error and is treated as 0.
797	 */
798	switch (up->dcycles & 0xff) {
799
800	case 0x00:		/* 0-1 bits (0) */
801	case 0x80:
802		irig_decode(peer, BIT0);
803		break;
804
805	case 0xc0:		/* 2-4 bits (1) */
806	case 0xe0:
807	case 0xf0:
808		irig_decode(peer, BIT1);
809		break;
810
811	case 0xf8:		/* (5-7 bits (PI) */
812	case 0xfc:
813	case 0xfe:
814		irig_decode(peer, BITP);
815		break;
816
817	default:		/* 8 bits (error) */
818		irig_decode(peer, BIT0);
819		up->errflg |= IRIG_ERR_DECODE;
820	}
821}
822
823
824/*
825 * irig_decode - decode the data
826 *
827 * This routine assembles bauds into digits, digits into frames and
828 * frames into the timecode fields. Bits can have values of zero, one
829 * or position identifier. There are four bits per digit, ten digits per
830 * frame and ten frames per second.
831 */
832static void
833irig_decode(
834	struct	peer *peer,	/* peer structure pointer */
835	int	bit		/* data bit (0, 1 or 2) */
836	)
837{
838	struct refclockproc *pp;
839	struct irigunit *up;
840
841	/*
842	 * Local variables
843	 */
844	int	syncdig;	/* sync digit (Spectracom) */
845	char	sbs[6 + 1];	/* binary seconds since 0h */
846	char	spare[2 + 1];	/* mulligan digits */
847	int	temp;
848
849	syncdig = 0;
850	pp = peer->procptr;
851	up = pp->unitptr;
852
853	/*
854	 * Assemble frame bits.
855	 */
856	up->bits >>= 1;
857	if (bit == BIT1) {
858		up->bits |= 0x200;
859	} else if (bit == BITP && up->lastbit == BITP) {
860
861		/*
862		 * Frame sync - two adjacent position identifiers, which
863		 * mark the beginning of the second. The reference time
864		 * is the beginning of the second position identifier,
865		 * so copy the character timestamp to the reference
866		 * timestamp.
867		 */
868		if (up->frmcnt != 1)
869			up->errflg |= IRIG_ERR_SYNCH;
870		up->frmcnt = 1;
871		up->refstamp = up->prvstamp;
872	}
873	up->lastbit = bit;
874	if (up->frmcnt % SUBFLD == 0) {
875
876		/*
877		 * End of frame. Encode two hexadecimal digits in
878		 * little-endian timecode field. Note frame 1 is shifted
879		 * right one bit to account for the marker PI.
880		 */
881		temp = up->bits;
882		if (up->frmcnt == 10)
883			temp >>= 1;
884		if (up->xptr >= 2) {
885			up->timecode[--up->xptr] = hexchar[temp & 0xf];
886			up->timecode[--up->xptr] = hexchar[(temp >> 5) &
887			    0xf];
888		}
889		if (up->frmcnt == 0) {
890
891			/*
892			 * End of second. Decode the timecode and wind
893			 * the clock. Not all IRIG generators have the
894			 * year; if so, it is nonzero after year 2000.
895			 * Not all have the hardware status bit; if so,
896			 * it is lit when the source is okay and dim
897			 * when bad. We watch this only if the year is
898			 * nonzero. Not all are configured for signature
899			 * control. If so, all BCD digits are set to
900			 * zero if the source is bad. In this case the
901			 * refclock_process() will reject the timecode
902			 * as invalid.
903			 */
904			up->xptr = 2 * SUBFLD;
905			if (sscanf((char *)up->timecode,
906			   "%6s%2d%1d%2s%3d%2d%2d%2d", sbs, &pp->year,
907			    &syncdig, spare, &pp->day, &pp->hour,
908			    &pp->minute, &pp->second) != 8)
909				pp->leap = LEAP_NOTINSYNC;
910			else
911				pp->leap = LEAP_NOWARNING;
912			up->second = (up->second + up->decim) % 60;
913
914			/*
915			 * Raise an alarm if the day field is zero,
916			 * which happens when signature control is
917			 * enabled and the device has lost
918			 * synchronization. Raise an alarm if the year
919			 * field is nonzero and the sync indicator is
920			 * zero, which happens when a Spectracom radio
921			 * has lost synchronization. Raise an alarm if
922			 * the expected second does not agree with the
923			 * decoded second, which happens with a garbled
924			 * IRIG signal. We are very particular.
925			 */
926			if (pp->day == 0 || (pp->year != 0 && syncdig ==
927			    0))
928				up->errflg |= IRIG_ERR_SIGERR;
929			if (pp->second != up->second)
930				up->errflg |= IRIG_ERR_CHECK;
931			up->second = pp->second;
932
933			/*
934			 * Wind the clock only if there are no errors
935			 * and the time constant has reached the
936			 * maximum.
937			 */
938			if (up->errflg == 0 && up->tc == MAXTC) {
939				pp->lastref = pp->lastrec;
940				pp->lastrec = up->refstamp;
941				if (!refclock_process(pp))
942					refclock_report(peer,
943					    CEVNT_BADTIME);
944			}
945			snprintf(pp->a_lastcode, sizeof(pp->a_lastcode),
946			    "%02x %02d %03d %02d:%02d:%02d %4.0f %3d %6.3f %2d %6.2f %6.1f %s",
947			    up->errflg, pp->year, pp->day,
948			    pp->hour, pp->minute, pp->second,
949			    up->maxsignal, up->gain, up->modndx,
950			    up->tc, up->exing * 1e6 / SECOND, up->freq *
951			    1e6 / SECOND, ulfptoa(&pp->lastrec, 6));
952			pp->lencode = strlen(pp->a_lastcode);
953			up->errflg = 0;
954			if (pp->sloppyclockflag & CLK_FLAG4) {
955				record_clock_stats(&peer->srcadr,
956				    pp->a_lastcode);
957#ifdef DEBUG
958				if (debug)
959					printf("irig %s\n",
960					    pp->a_lastcode);
961#endif /* DEBUG */
962			}
963		}
964	}
965	up->frmcnt = (up->frmcnt + 1) % FIELD;
966}
967
968
969/*
970 * irig_poll - called by the transmit procedure
971 *
972 * This routine sweeps up the timecode updates since the last poll. For
973 * IRIG-B there should be at least 60 updates; for IRIG-E there should
974 * be at least 6. If nothing is heard, a timeout event is declared.
975 */
976static void
977irig_poll(
978	int	unit,		/* instance number (not used) */
979	struct peer *peer	/* peer structure pointer */
980	)
981{
982	struct refclockproc *pp;
983
984	pp = peer->procptr;
985
986	if (pp->coderecv == pp->codeproc) {
987		refclock_report(peer, CEVNT_TIMEOUT);
988		return;
989
990	}
991	refclock_receive(peer);
992	if (!(pp->sloppyclockflag & CLK_FLAG4)) {
993		record_clock_stats(&peer->srcadr, pp->a_lastcode);
994#ifdef DEBUG
995		if (debug)
996			printf("irig %s\n", pp->a_lastcode);
997#endif /* DEBUG */
998	}
999	pp->polls++;
1000
1001}
1002
1003
1004/*
1005 * irig_gain - adjust codec gain
1006 *
1007 * This routine is called at the end of each second. It uses the AGC to
1008 * bradket the maximum signal level between MINAMP and MAXAMP to avoid
1009 * hunting. The routine also jiggles the input port and selectively
1010 * mutes the monitor.
1011 */
1012static void
1013irig_gain(
1014	struct peer *peer	/* peer structure pointer */
1015	)
1016{
1017	struct refclockproc *pp;
1018	struct irigunit *up;
1019
1020	pp = peer->procptr;
1021	up = pp->unitptr;
1022
1023	/*
1024	 * Apparently, the codec uses only the high order bits of the
1025	 * gain control field. Thus, it may take awhile for changes to
1026	 * wiggle the hardware bits.
1027	 */
1028	if (up->maxsignal < MINAMP) {
1029		up->gain += 4;
1030		if (up->gain > MAXGAIN)
1031			up->gain = MAXGAIN;
1032	} else if (up->maxsignal > MAXAMP) {
1033		up->gain -= 4;
1034		if (up->gain < 0)
1035			up->gain = 0;
1036	}
1037	audio_gain(up->gain, up->mongain, up->port);
1038}
1039
1040
1041#else
1042int refclock_irig_bs;
1043#endif /* REFCLOCK */
1044