/*
* Audio WWV/H demodulator/decoder
*
* This driver synchronizes the computer time using data encoded in
* radio transmissions from NIST time/frequency stations WWV in Boulder,
* CO, and WWVH in Kauai, HI. Transmissions are made continuously on
* 2.5, 5, 10 and 15 MHz from WWV and WWVH, and 20 MHz from WWV. An
* ordinary AM shortwave receiver can be tuned manually to one of these
* frequencies or, in the case of ICOM receivers, the receiver can be
* tuned automatically using this program as propagation conditions
* change throughout the weasons, both day and night.
*
* The driver requires an audio codec or sound card with sampling rate 8
* kHz and mu-law companding. This is the same standard as used by the
* telephone industry and is supported by most hardware and operating
* systems, including Solaris, SunOS, FreeBSD, NetBSD and Linux. In this
* implementation, only one audio driver and codec can be supported on a
* single machine.
*
* The demodulation and decoding algorithms used in this driver are
* based on those developed for the TAPR DSP93 development board and the
* TI 320C25 digital signal processor described in: Mills, D.L. A
* precision radio clock for WWV transmissions. Electrical Engineering
* Report 97-8-1, University of Delaware, August 1997, 25 pp., available
* from www.eecis.udel.edu/~mills/reports.html. The algorithms described
* in this report have been modified somewhat to improve performance
* under weak signal conditions and to provide an automatic station
* identification feature.
*
* The ICOM code is normally compiled in the driver. It isn't used,
* unless the mode keyword on the server configuration command specifies
* a nonzero ICOM ID select code. The C-IV trace is turned on if the
* debug level is greater than one.
*
* Fudge factors
*
* Fudge flag4 causes the debugging output described above to be
* recorded in the clockstats file. Fudge flag2 selects the audio input
* port, where 0 is the mike port (default) and 1 is the line-in port.
* It does not seem useful to select the compact disc player port. Fudge
* flag3 enables audio monitoring of the input signal. For this purpose,
* the monitor gain is set to a default value.
*
* CEVNT_BADTIME invalid date or time
* CEVNT_PROP propagation failure - no stations heard
* CEVNT_TIMEOUT timeout (see newgame() below)
*/
/*
* General definitions. These ordinarily do not need to be changed.
*/
#define DEVICE_AUDIO "/dev/audio" /* audio device name */
#define AUDIO_BUFSIZ 320 /* audio buffer size (50 ms) */
#define PRECISION (-10) /* precision assumed (about 1 ms) */
#define DESCRIPTION "WWV/H Audio Demodulator/Decoder" /* WRU */
#define WWV_SEC 8000 /* second epoch (sample rate) (Hz) */
#define WWV_MIN (WWV_SEC * 60) /* minute epoch */
#define OFFSET 128 /* companded sample offset */
#define SIZE 256 /* decompanding table size */
#define MAXAMP 6000. /* max signal level reference */
#define MAXCLP 100 /* max clips above reference per s */
#define MAXSNR 40. /* max SNR reference */
#define MAXFREQ 1.5 /* max frequency tolerance (187 PPM) */
#define DATCYC 170 /* data filter cycles */
#define DATSIZ (DATCYC * MS) /* data filter size */
#define SYNCYC 800 /* minute filter cycles */
#define SYNSIZ (SYNCYC * MS) /* minute filter size */
#define TCKCYC 5 /* tick filter cycles */
#define TCKSIZ (TCKCYC * MS) /* tick filter size */
#define NCHAN 5 /* number of radio channels */
#define AUDIO_PHI 5e-6 /* dispersion growth factor */
#define TBUF 128 /* max monitor line length */
/*
* Tunable parameters. The DGAIN parameter can be changed to fit the
* audio response of the radio at 100 Hz. The WWV/WWVH data subcarrier
* is transmitted at about 20 percent percent modulation; the matched
* filter boosts it by a factor of 17 and the receiver response does
* what it does. The compromise value works for ICOM radios. If the
* radio is not tunable, the DCHAN parameter can be changed to fit the
* expected best propagation frequency: higher if further from the
* transmitter, lower if nearer. The compromise value works for the US
* right coast.
*/
#define DCHAN 3 /* default radio channel (15 Mhz) */
#define DGAIN 5. /* subcarrier gain */
/*
* General purpose status bits (status)
*
* SELV and/or SELH are set when WWV or WWVH have been heard and cleared
* on signal loss. SSYNC is set when the second sync pulse has been
* acquired and cleared by signal loss. MSYNC is set when the minute
* sync pulse has been acquired. DSYNC is set when the units digit has
* has reached the threshold and INSYNC is set when all nine digits have
* reached the threshold. The MSYNC, DSYNC and INSYNC bits are cleared
* only by timeout, upon which the driver starts over from scratch.
*
* DGATE is lit if the data bit amplitude or SNR is below thresholds and
* BGATE is lit if the pulse width amplitude or SNR is below thresolds.
* LEPSEC is set during the last minute of the leap day. At the end of
* this minute the driver inserts second 60 in the seconds state machine
* and the minute sync slips a second.
*/
#define MSYNC 0x0001 /* minute epoch sync */
#define SSYNC 0x0002 /* second epoch sync */
#define DSYNC 0x0004 /* minute units sync */
#define INSYNC 0x0008 /* clock synchronized */
#define FGATE 0x0010 /* frequency gate */
#define DGATE 0x0020 /* data pulse amplitude error */
#define BGATE 0x0040 /* data pulse width error */
#define METRIC 0x0080 /* one or more stations heard */
#define LEPSEC 0x1000 /* leap minute */
/*
* Station scoreboard bits
*
* These are used to establish the signal quality for each of the five
* frequencies and two stations.
*/
#define SELV 0x0100 /* WWV station select */
#define SELH 0x0200 /* WWVH station select */
/*
* Alarm status bits (alarm)
*
* These bits indicate various alarm conditions, which are decoded to
* form the quality character included in the timecode.
*/
#define CMPERR 0x1 /* digit or misc bit compare error */
#define LOWERR 0x2 /* low bit or digit amplitude or SNR */
#define NINERR 0x4 /* less than nine digits in minute */
#define SYNERR 0x8 /* not tracking second sync */
/*
* Watchcat timeouts (watch)
*
* If these timeouts expire, the status bits are mashed to zero and the
* driver starts from scratch. Suitably more refined procedures may be
* developed in future. All these are in minutes.
*/
#define ACQSN 6 /* station acquisition timeout */
#define DATA 15 /* unit minutes timeout */
#define SYNCH 40 /* station sync timeout */
#define PANIC (2 * 1440) /* panic timeout */
/*
* Thresholds. These establish the minimum signal level, minimum SNR and
* maximum jitter thresholds which establish the error and false alarm
* rates of the driver. The values defined here may be on the
* adventurous side in the interest of the highest sensitivity.
*/
#define MTHR 13. /* minute sync gate (percent) */
#define TTHR 50. /* minute sync threshold (percent) */
#define AWND 20 /* minute sync jitter threshold (ms) */
#define ATHR 2500. /* QRZ minute sync threshold */
#define ASNR 20. /* QRZ minute sync SNR threshold (dB) */
#define QTHR 2500. /* QSY minute sync threshold */
#define QSNR 20. /* QSY minute sync SNR threshold (dB) */
#define STHR 2500. /* second sync threshold */
#define SSNR 15. /* second sync SNR threshold (dB) */
#define SCMP 10 /* second sync compare threshold */
#define DTHR 1000. /* bit threshold */
#define DSNR 10. /* bit SNR threshold (dB) */
#define AMIN 3 /* min bit count */
#define AMAX 6 /* max bit count */
#define BTHR 1000. /* digit threshold */
#define BSNR 3. /* digit likelihood threshold (dB) */
#define BCMP 3 /* digit compare threshold */
#define MAXERR 40 /* maximum error alarm */
/*
* Tone frequency definitions. The increments are for 4.5-deg sine
* table.
*/
#define MS (WWV_SEC / 1000) /* samples per millisecond */
#define IN100 ((100 * 80) / WWV_SEC) /* 100 Hz increment */
#define IN1000 ((1000 * 80) / WWV_SEC) /* 1000 Hz increment */
#define IN1200 ((1200 * 80) / WWV_SEC) /* 1200 Hz increment */
/*
* Acquisition and tracking time constants
*/
#define MINAVG 8 /* min averaging time */
#define MAXAVG 1024 /* max averaging time */
#define FCONST 3 /* frequency time constant */
#define TCONST 16 /* data bit/digit time constant */
/*
* Miscellaneous status bits (misc)
*
* These bits correspond to designated bits in the WWV/H timecode. The
* bit probabilities are exponentially averaged over several minutes and
* processed by a integrator and threshold.
*/
#define DUT1 0x01 /* 56 DUT .1 */
#define DUT2 0x02 /* 57 DUT .2 */
#define DUT4 0x04 /* 58 DUT .4 */
#define DUTS 0x08 /* 50 DUT sign */
#define DST1 0x10 /* 55 DST1 leap warning */
#define DST2 0x20 /* 2 DST2 DST1 delayed one day */
#define SECWAR 0x40 /* 3 leap second warning */
/*
* The on-time synchronization point is the positive-going zero crossing
* of the first cycle of the 5-ms second pulse. The IIR baseband filter
* phase delay is 0.91 ms, while the receiver delay is approximately 4.7
* ms at 1000 Hz. The fudge value -0.45 ms due to the codec and other
* causes was determined by calibrating to a PPS signal from a GPS
* receiver. The additional propagation delay specific to each receiver
* location can be programmed in the fudge time1 and time2 values for
* WWV and WWVH, respectively.
*
* The resulting offsets with a 2.4-GHz P4 running FreeBSD 6.1 are
* generally within .02 ms short-term with .02 ms jitter. The long-term
* offsets vary up to 0.3 ms due to ionosperhic layer height variations.
* The processor load due to the driver is 5.8 percent.
*/
#define PDELAY ((.91 + 4.7 - 0.45) / 1000) /* system delay (s) */
/*
* Decoder operations at the end of each second are driven by a state
* machine. The transition matrix consists of a dispatch table indexed
* by second number. Each entry in the table contains a case switch
* number and argument.
*/
struct progx {
int sw; /* case switch number */
int arg; /* argument */
};
/*
* Case switch numbers
*/
#define IDLE 0 /* no operation */
#define COEF 1 /* BCD bit */
#define COEF1 2 /* BCD bit for minute unit */
#define COEF2 3 /* BCD bit not used */
#define DECIM9 4 /* BCD digit 0-9 */
#define DECIM6 5 /* BCD digit 0-6 */
#define DECIM3 6 /* BCD digit 0-3 */
#define DECIM2 7 /* BCD digit 0-2 */
#define MSCBIT 8 /* miscellaneous bit */
#define MSC20 9 /* miscellaneous bit */
#define MSC21 10 /* QSY probe channel */
#define MIN1 11 /* latch time */
#define MIN2 12 /* leap second */
#define SYNC2 13 /* latch minute sync pulse */
#define SYNC3 14 /* latch data pulse */
/*
* Offsets in decoding matrix
*/
#define MN 0 /* minute digits (2) */
#define HR 2 /* hour digits (2) */
#define DA 4 /* day digits (3) */
#define YR 7 /* year digits (2) */
/*
* DST decode (DST2 DST1) for prettyprint
*/
char dstcod[] = {
'S', /* 00 standard time */
'I', /* 01 set clock ahead at 0200 local */
'O', /* 10 set clock back at 0200 local */
'D' /* 11 daylight time */
};
/*
* The decoding matrix consists of nine row vectors, one for each digit
* of the timecode. The digits are stored from least to most significant
* order. The maximum-likelihood timecode is formed from the digits
* corresponding to the maximum-likelihood values reading in the
* opposite order: yy ddd hh:mm.
*/
struct decvec {
int radix; /* radix (3, 4, 6, 10) */
int digit; /* current clock digit */
int count; /* match count */
double digprb; /* max digit probability */
double digsnr; /* likelihood function (dB) */
double like[10]; /* likelihood integrator 0-9 */
};
/*
* The station structure (sp) is used to acquire the minute pulse from
* WWV and/or WWVH. These stations are distinguished by the frequency
* used for the second and minute sync pulses, 1000 Hz for WWV and 1200
* Hz for WWVH. Other than frequency, the format is the same.
*/
struct sync {
double epoch; /* accumulated epoch differences */
double maxeng; /* sync max energy */
double noieng; /* sync noise energy */
long pos; /* max amplitude position */
long lastpos; /* last max position */
long mepoch; /* minute synch epoch */
double amp; /* sync signal */
double syneng; /* sync signal max */
double synmax; /* sync signal max latched at 0 s */
double synsnr; /* sync signal SNR */
double metric; /* signal quality metric */
int reach; /* reachability register */
int count; /* bit counter */
int select; /* select bits */
char refid[5]; /* reference identifier */
};
/*
* The channel structure (cp) is used to mitigate between channels.
*/
struct chan {
int gain; /* audio gain */
struct sync wwv; /* wwv station */
struct sync wwvh; /* wwvh station */
};
/*
* WWV unit control structure (up)
*/
struct wwvunit {
l_fp timestamp; /* audio sample timestamp */
l_fp tick; /* audio sample increment */
double phase, freq; /* logical clock phase and frequency */
double monitor; /* audio monitor point */
double pdelay; /* propagation delay (s) */
#ifdef ICOM
int fd_icom; /* ICOM file descriptor */
#endif /* ICOM */
int errflg; /* error flags */
int watch; /* watchcat */
/*
* Audio codec variables
*/
double comp[SIZE]; /* decompanding table */
int port; /* codec port */
int gain; /* codec gain */
int mongain; /* codec monitor gain */
int clipcnt; /* sample clipped count */
/*
* Variables used to establish basic system timing
*/
int avgint; /* master time constant */
int yepoch; /* sync epoch */
int repoch; /* buffered sync epoch */
double epomax; /* second sync amplitude */
double eposnr; /* second sync SNR */
double irig; /* data I channel amplitude */
double qrig; /* data Q channel amplitude */
int datapt; /* 100 Hz ramp */
double datpha; /* 100 Hz VFO control */
int rphase; /* second sample counter */
long mphase; /* minute sample counter */
/*
* Variables used to mitigate which channel to use
*/
struct chan mitig[NCHAN]; /* channel data */
struct sync *sptr; /* station pointer */
int dchan; /* data channel */
int schan; /* probe channel */
int achan; /* active channel */
/*
* Variables used by the clock state machine
*/
struct decvec decvec[9]; /* decoding matrix */
int rsec; /* seconds counter */
int digcnt; /* count of digits synchronized */
/*
* Variables used to estimate signal levels and bit/digit
* probabilities
*/
double datsig; /* data signal max */
double datsnr; /* data signal SNR (dB) */
/*
* Variables used to establish status and alarm conditions
*/
int status; /* status bits */
int alarm; /* alarm flashers */
int misc; /* miscellaneous timecode bits */
int errcnt; /* data bit error counter */
};
/*
* Transfer vector
*/
struct refclock refclock_wwv = {
wwv_start, /* start up driver */
wwv_shutdown, /* shut down driver */
wwv_poll, /* transmit poll message */
noentry, /* not used (old wwv_control) */
noentry, /* initialize driver (not used) */
noentry, /* not used (old wwv_buginfo) */
NOFLAGS /* not used */
};
/*
* wwv_start - open the devices and initialize data for processing
*/
static int
wwv_start(
int unit, /* instance number (used by PCM) */
struct peer *peer /* peer structure pointer */
)
{
struct refclockproc *pp;
struct wwvunit *up;
#ifdef ICOM
int temp;
#endif /* ICOM */
/*
* Local variables
*/
int fd; /* file descriptor */
int i; /* index */
double step; /* codec adjustment */
/*
* Open audio device
*/
fd = audio_init(DEVICE_AUDIO, AUDIO_BUFSIZ, unit);
if (fd < 0)
return (0);
#ifdef DEBUG
if (debug)
audio_show();
#endif /* DEBUG */
/*
* Allocate and initialize unit structure
*/
up = emalloc_zero(sizeof(*up));
pp = peer->procptr;
pp->io.clock_recv = wwv_receive;
pp->io.srcclock = peer;
pp->io.datalen = 0;
pp->io.fd = fd;
if (!io_addclock(&pp->io)) {
close(fd);
free(up);
return (0);
}
pp->unitptr = up;
/*
* The companded samples are encoded sign-magnitude. The table
* contains all the 256 values in the interest of speed.
*/
up->comp[0] = up->comp[OFFSET] = 0.;
up->comp[1] = 1.; up->comp[OFFSET + 1] = -1.;
up->comp[2] = 3.; up->comp[OFFSET + 2] = -3.;
step = 2.;
for (i = 3; i < OFFSET; i++) {
up->comp[i] = up->comp[i - 1] + step;
up->comp[OFFSET + i] = -up->comp[i];
if (i % 16 == 0)
step *= 2.;
}
DTOLFP(1. / WWV_SEC, &up->tick);
/*
* Initialize the decoding matrix with the radix for each digit
* position.
*/
up->decvec[MN].radix = 10; /* minutes */
up->decvec[MN + 1].radix = 6;
up->decvec[HR].radix = 10; /* hours */
up->decvec[HR + 1].radix = 3;
up->decvec[DA].radix = 10; /* days */
up->decvec[DA + 1].radix = 10;
up->decvec[DA + 2].radix = 4;
up->decvec[YR].radix = 10; /* years */
up->decvec[YR + 1].radix = 10;
#ifdef ICOM
/*
* Initialize autotune if available. Note that the ICOM select
* code must be less than 128, so the high order bit can be used
* to select the line speed 0 (9600 bps) or 1 (1200 bps). Note
* we don't complain if the ICOM device is not there; but, if it
* is, the radio better be working.
*/
temp = 0;
#ifdef DEBUG
if (debug > 1)
temp = P_TRACE;
#endif /* DEBUG */
if (peer->ttl != 0) {
if (peer->ttl & 0x80)
up->fd_icom = icom_init("/dev/icom", B1200,
temp);
else
up->fd_icom = icom_init("/dev/icom", B9600,
temp);
}
if (up->fd_icom > 0) {
if (wwv_qsy(peer, DCHAN) != 0) {
msyslog(LOG_NOTICE, "icom: radio not found");
close(up->fd_icom);
up->fd_icom = 0;
} else {
msyslog(LOG_NOTICE, "icom: autotune enabled");
}
}
#endif /* ICOM */
/*
* Let the games begin.
*/
wwv_newgame(peer);
return (1);
}
/*
* wwv_shutdown - shut down the clock
*/
static void
wwv_shutdown(
int unit, /* instance number (not used) */
struct peer *peer /* peer structure pointer */
)
{
struct refclockproc *pp;
struct wwvunit *up;
pp = peer->procptr;
up = pp->unitptr;
if (up == NULL)
return;
/*
* wwv_receive - receive data from the audio device
*
* This routine reads input samples and adjusts the logical clock to
* track the A/D sample clock by dropping or duplicating codec samples.
* It also controls the A/D signal level with an AGC loop to mimimize
* quantization noise and avoid overload.
*/
static void
wwv_receive(
struct recvbuf *rbufp /* receive buffer structure pointer */
)
{
struct peer *peer;
struct refclockproc *pp;
struct wwvunit *up;
peer = rbufp->recv_peer;
pp = peer->procptr;
up = pp->unitptr;
/*
* Main loop - read until there ain't no more. Note codec
* samples are bit-inverted.
*/
DTOLFP((double)rbufp->recv_length / WWV_SEC, <emp);
L_SUB(&rbufp->recv_time, <emp);
up->timestamp = rbufp->recv_time;
dpt = rbufp->recv_buffer;
for (bufcnt = 0; bufcnt < rbufp->recv_length; bufcnt++) {
sample = up->comp[~*dpt++ & 0xff];
/*
* Clip noise spikes greater than MAXAMP (6000) and
* record the number of clips to be used later by the
* AGC.
*/
if (sample > MAXAMP) {
sample = MAXAMP;
up->clipcnt++;
} else if (sample < -MAXAMP) {
sample = -MAXAMP;
up->clipcnt++;
}
/*
* Variable frequency oscillator. The codec oscillator
* runs at the nominal rate of 8000 samples per second,
* or 125 us per sample. A frequency change of one unit
* results in either duplicating or deleting one sample
* per second, which results in a frequency change of
* 125 PPM.
*/
up->phase += (up->freq + clock_codec) / WWV_SEC;
if (up->phase >= .5) {
up->phase -= 1.;
} else if (up->phase < -.5) {
up->phase += 1.;
wwv_rf(peer, sample);
wwv_rf(peer, sample);
} else {
wwv_rf(peer, sample);
}
L_ADD(&up->timestamp, &up->tick);
}
/*
* Set the input port and monitor gain for the next buffer.
*/
if (pp->sloppyclockflag & CLK_FLAG2)
up->port = 2;
else
up->port = 1;
if (pp->sloppyclockflag & CLK_FLAG3)
up->mongain = MONGAIN;
else
up->mongain = 0;
}
/*
* wwv_poll - called by the transmit procedure
*
* This routine keeps track of status. If no offset samples have been
* processed during a poll interval, a timeout event is declared. If
* errors have have occurred during the interval, they are reported as
* well.
*/
static void
wwv_poll(
int unit, /* instance number (not used) */
struct peer *peer /* peer structure pointer */
)
{
struct refclockproc *pp;
struct wwvunit *up;
pp = peer->procptr;
up = pp->unitptr;
if (up->errflg)
refclock_report(peer, up->errflg);
up->errflg = 0;
pp->polls++;
}
/*
* wwv_rf - process signals and demodulate to baseband
*
* This routine grooms and filters decompanded raw audio samples. The
* output signal is the 100-Hz filtered baseband data signal in
* quadrature phase. The routine also determines the minute synch epoch,
* as well as certain signal maxima, minima and related values.
*
* There are two 1-s ramps used by this program. Both count the 8000
* logical clock samples spanning exactly one second. The epoch ramp
* counts the samples starting at an arbitrary time. The rphase ramp
* counts the samples starting at the 5-ms second sync pulse found
* during the epoch ramp.
*
* There are two 1-m ramps used by this program. The mphase ramp counts
* the 480,000 logical clock samples spanning exactly one minute and
* starting at an arbitrary time. The rsec ramp counts the 60 seconds of
* the minute starting at the 800-ms minute sync pulse found during the
* mphase ramp. The rsec ramp drives the seconds state machine to
* determine the bits and digits of the timecode.
*
* Demodulation operations are based on three synthesized quadrature
* sinusoids: 100 Hz for the data signal, 1000 Hz for the WWV sync
* signal and 1200 Hz for the WWVH sync signal. These drive synchronous
* matched filters for the data signal (170 ms at 100 Hz), WWV minute
* sync signal (800 ms at 1000 Hz) and WWVH minute sync signal (800 ms
* at 1200 Hz). Two additional matched filters are switched in
* as required for the WWV second sync signal (5 cycles at 1000 Hz) and
* WWVH second sync signal (6 cycles at 1200 Hz).
*/
static void
wwv_rf(
struct peer *peer, /* peerstructure pointer */
double isig /* input signal */
)
{
struct refclockproc *pp;
struct wwvunit *up;
struct sync *sp, *rp;
static int iptr; /* data channel pointer */
static double ibuf[DATSIZ]; /* data I channel delay line */
static double qbuf[DATSIZ]; /* data Q channel delay line */
static int jptr; /* sync channel pointer */
static int kptr; /* tick channel pointer */
static int csinptr; /* wwv channel phase */
static double cibuf[SYNSIZ]; /* wwv I channel delay line */
static double cqbuf[SYNSIZ]; /* wwv Q channel delay line */
static double ciamp; /* wwv I channel amplitude */
static double cqamp; /* wwv Q channel amplitude */
static double csibuf[TCKSIZ]; /* wwv I tick delay line */
static double csqbuf[TCKSIZ]; /* wwv Q tick delay line */
static double csiamp; /* wwv I tick amplitude */
static double csqamp; /* wwv Q tick amplitude */
static int hsinptr; /* wwvh channel phase */
static double hibuf[SYNSIZ]; /* wwvh I channel delay line */
static double hqbuf[SYNSIZ]; /* wwvh Q channel delay line */
static double hiamp; /* wwvh I channel amplitude */
static double hqamp; /* wwvh Q channel amplitude */
static double hsibuf[TCKSIZ]; /* wwvh I tick delay line */
static double hsqbuf[TCKSIZ]; /* wwvh Q tick delay line */
static double hsiamp; /* wwvh I tick amplitude */
static double hsqamp; /* wwvh Q tick amplitude */
static double epobuf[WWV_SEC]; /* second sync comb filter */
static double epomax, nxtmax; /* second sync amplitude buffer */
static int epopos; /* epoch second sync position buffer */
static int iniflg; /* initialization flag */
int epoch; /* comb filter index */
double dtemp;
int i;
/*
* Baseband data demodulation. The 100-Hz subcarrier is
* extracted using a 150-Hz IIR lowpass filter. This attenuates
* the 1000/1200-Hz sync signals, as well as the 440-Hz and
* 600-Hz tones and most of the noise and voice modulation
* components.
*
* The subcarrier is transmitted 10 dB down from the carrier.
* The DGAIN parameter can be adjusted for this and to
* compensate for the radio audio response at 100 Hz.
*
* Matlab IIR 4th-order IIR elliptic, 150 Hz lowpass, 0.2 dB
* passband ripple, -50 dB stopband ripple, phase delay 0.97 ms.
*/
data = (lpf[4] = lpf[3]) * 8.360961e-01;
data += (lpf[3] = lpf[2]) * -3.481740e+00;
data += (lpf[2] = lpf[1]) * 5.452988e+00;
data += (lpf[1] = lpf[0]) * -3.807229e+00;
lpf[0] = isig * DGAIN - data;
data = lpf[0] * 3.281435e-03
+ lpf[1] * -1.149947e-02
+ lpf[2] * 1.654858e-02
+ lpf[3] * -1.149947e-02
+ lpf[4] * 3.281435e-03;
/*
* The 100-Hz data signal is demodulated using a pair of
* quadrature multipliers, matched filters and a phase lock
* loop. The I and Q quadrature data signals are produced by
* multiplying the filtered signal by 100-Hz sine and cosine
* signals, respectively. The signals are processed by 170-ms
* synchronous matched filters to produce the amplitude and
* phase signals used by the demodulator. The signals are scaled
* to produce unit energy at the maximum value.
*/
i = up->datapt;
up->datapt = (up->datapt + IN100) % 80;
dtemp = sintab[i] * data / (MS / 2. * DATCYC);
up->irig -= ibuf[iptr];
ibuf[iptr] = dtemp;
up->irig += dtemp;
/*
* Baseband sync demodulation. The 1000/1200 sync signals are
* extracted using a 600-Hz IIR bandpass filter. This removes
* the 100-Hz data subcarrier, as well as the 440-Hz and 600-Hz
* tones and most of the noise and voice modulation components.
*
* Matlab 4th-order IIR elliptic, 800-1400 Hz bandpass, 0.2 dB
* passband ripple, -50 dB stopband ripple, phase delay 0.91 ms.
*/
syncx = (bpf[8] = bpf[7]) * 4.897278e-01;
syncx += (bpf[7] = bpf[6]) * -2.765914e+00;
syncx += (bpf[6] = bpf[5]) * 8.110921e+00;
syncx += (bpf[5] = bpf[4]) * -1.517732e+01;
syncx += (bpf[4] = bpf[3]) * 1.975197e+01;
syncx += (bpf[3] = bpf[2]) * -1.814365e+01;
syncx += (bpf[2] = bpf[1]) * 1.159783e+01;
syncx += (bpf[1] = bpf[0]) * -4.735040e+00;
bpf[0] = isig - syncx;
syncx = bpf[0] * 8.203628e-03
+ bpf[1] * -2.375732e-02
+ bpf[2] * 3.353214e-02
+ bpf[3] * -4.080258e-02
+ bpf[4] * 4.605479e-02
+ bpf[5] * -4.080258e-02
+ bpf[6] * 3.353214e-02
+ bpf[7] * -2.375732e-02
+ bpf[8] * 8.203628e-03;
/*
* The 1000/1200 sync signals are demodulated using a pair of
* quadrature multipliers and matched filters. However,
* synchronous demodulation at these frequencies is impractical,
* so only the signal amplitude is used. The I and Q quadrature
* sync signals are produced by multiplying the filtered signal
* by 1000-Hz (WWV) and 1200-Hz (WWVH) sine and cosine signals,
* respectively. The WWV and WWVH signals are processed by 800-
* ms synchronous matched filters and combined to produce the
* minute sync signal and detect which one (or both) the WWV or
* WWVH signal is present. The WWV and WWVH signals are also
* processed by 5-ms synchronous matched filters and combined to
* produce the second sync signal. The signals are scaled to
* produce unit energy at the maximum value.
*
* Note the master timing ramps, which run continuously. The
* minute counter (mphase) counts the samples in the minute,
* while the second counter (epoch) counts the samples in the
* second.
*/
up->mphase = (up->mphase + 1) % WWV_MIN;
epoch = up->mphase % WWV_SEC;
/*
* The following section is called once per minute. It does
* housekeeping and timeout functions and empties the dustbins.
*/
if (up->mphase == 0) {
up->watch++;
if (!(up->status & MSYNC)) {
/*
* If minute sync has not been acquired before
* ACQSN timeout (6 min), or if no signal is
* heard, the program cycles to the next
* frequency and tries again.
*/
if (!wwv_newchan(peer))
up->watch = 0;
} else {
/*
* If the leap bit is set, set the minute epoch
* back one second so the station processes
* don't miss a beat.
*/
if (up->status & LEPSEC) {
up->mphase -= WWV_SEC;
if (up->mphase < 0)
up->mphase += WWV_MIN;
}
}
}
/*
* When the channel metric reaches threshold and the second
* counter matches the minute epoch within the second, the
* driver has synchronized to the station. The second number is
* the remaining seconds until the next minute epoch, while the
* sync epoch is zero. Watch out for the first second; if
* already synchronized to the second, the buffered sync epoch
* must be set.
*
* Note the guard interval is 200 ms; if for some reason the
* clock drifts more than that, it might wind up in the wrong
* second. If the maximum frequency error is not more than about
* 1 PPM, the clock can go as much as two days while still in
* the same second.
*/
if (up->status & MSYNC) {
wwv_epoch(peer);
} else if (up->sptr != NULL) {
sp = up->sptr;
if (sp->metric >= TTHR && epoch == sp->mepoch % WWV_SEC)
{
up->rsec = (60 - sp->mepoch / WWV_SEC) % 60;
up->rphase = 0;
up->status |= MSYNC;
up->watch = 0;
if (!(up->status & SSYNC))
up->repoch = up->yepoch = epoch;
else
up->repoch = up->yepoch;
}
}
/*
* The second sync pulse is extracted using 5-ms (40 sample) FIR
* matched filters at 1000 Hz for WWV or 1200 Hz for WWVH. This
* pulse is used for the most precise synchronization, since if
* provides a resolution of one sample (125 us). The filters run
* only if the station has been reliably determined.
*/
if (up->status & SELV)
mfsync = sqrt(csiamp * csiamp + csqamp * csqamp) /
TCKCYC;
else if (up->status & SELH)
mfsync = sqrt(hsiamp * hsiamp + hsqamp * hsqamp) /
TCKCYC;
else
mfsync = 0;
/*
* Enhance the seconds sync pulse using a 1-s (8000-sample) comb
* filter. Correct for the FIR matched filter delay, which is 5
* ms for both the WWV and WWVH filters, and also for the
* propagation delay. Once each second look for second sync. If
* not in minute sync, fiddle the codec gain. Note the SNR is
* computed from the maximum sample and the envelope of the
* sample 6 ms before it, so if we slip more than a cycle the
* SNR should plummet. The signal is scaled to produce unit
* energy at the maximum value.
*/
dtemp = (epobuf[epoch] += (mfsync - epobuf[epoch]) /
up->avgint);
if (dtemp > epomax) {
int j;
/*
* wwv_qrz - identify and acquire WWV/WWVH minute sync pulse
*
* This routine implements a virtual station process used to acquire
* minute sync and to mitigate among the ten frequency and station
* combinations. During minute sync acquisition the process probes each
* frequency and station in turn for the minute pulse, which
* involves searching through the entire 480,000-sample minute. The
* process finds the maximum signal and RMS noise plus signal. Then, the
* actual noise is determined by subtracting the energy of the matched
* filter.
*
* Students of radar receiver technology will discover this algorithm
* amounts to a range-gate discriminator. A valid pulse must have peak
* amplitude at least QTHR (2500) and SNR at least QSNR (20) dB and the
* difference between the current and previous epoch must be less than
* AWND (20 ms). Note that the discriminator peak occurs about 800 ms
* into the second, so the timing is retarded to the previous second
* epoch.
*/
static void
wwv_qrz(
struct peer *peer, /* peer structure pointer */
struct sync *sp, /* sync channel structure */
int pdelay /* propagation delay (samples) */
)
{
struct refclockproc *pp;
struct wwvunit *up;
char tbuf[TBUF]; /* monitor buffer */
long epoch;
pp = peer->procptr;
up = pp->unitptr;
/*
* Find the sample with peak amplitude, which defines the minute
* epoch. Accumulate all samples to determine the total noise
* energy.
*/
epoch = up->mphase - pdelay - SYNSIZ;
if (epoch < 0)
epoch += WWV_MIN;
if (sp->amp > sp->maxeng) {
sp->maxeng = sp->amp;
sp->pos = epoch;
}
sp->noieng += sp->amp;
/*
* At the end of the minute, determine the epoch of the minute
* sync pulse, as well as the difference between the current and
* previous epoches due to the intrinsic frequency error plus
* jitter. When calculating the SNR, subtract the pulse energy
* from the total noise energy and then normalize.
*/
if (up->mphase == 0) {
sp->synmax = sp->maxeng;
sp->synsnr = wwv_snr(sp->synmax, (sp->noieng -
sp->synmax) / WWV_MIN);
if (sp->count == 0)
sp->lastpos = sp->pos;
epoch = (sp->pos - sp->lastpos) % WWV_MIN;
sp->reach <<= 1;
if (sp->reach & (1 << AMAX))
sp->count--;
if (sp->synmax > ATHR && sp->synsnr > ASNR) {
if (labs(epoch) < AWND * MS) {
sp->reach |= 1;
sp->count++;
sp->mepoch = sp->lastpos = sp->pos;
} else if (sp->count == 1) {
sp->lastpos = sp->pos;
}
}
if (up->watch > ACQSN)
sp->metric = 0;
else
sp->metric = wwv_metric(sp);
if (pp->sloppyclockflag & CLK_FLAG4) {
snprintf(tbuf, sizeof(tbuf),
"wwv8 %04x %3d %s %04x %.0f %.0f/%.1f %ld %ld",
up->status, up->gain, sp->refid,
sp->reach & 0xffff, sp->metric, sp->synmax,
sp->synsnr, sp->pos % WWV_SEC, epoch);
record_clock_stats(&peer->srcadr, tbuf);
#ifdef DEBUG
if (debug)
printf("%s\n", tbuf);
#endif /* DEBUG */
}
sp->maxeng = sp->noieng = 0;
}
}
/*
* wwv_endpoc - identify and acquire second sync pulse
*
* This routine is called at the end of the second sync interval. It
* determines the second sync epoch position within the second and
* disciplines the sample clock using a frequency-lock loop (FLL).
*
* Second sync is determined in the RF input routine as the maximum
* over all 8000 samples in the second comb filter. To assure accurate
* and reliable time and frequency discipline, this routine performs a
* great deal of heavy-handed heuristic data filtering and grooming.
*/
static void
wwv_endpoc(
struct peer *peer, /* peer structure pointer */
int epopos /* epoch max position */
)
{
struct refclockproc *pp;
struct wwvunit *up;
static int epoch_mf[3]; /* epoch median filter */
static int tepoch; /* current second epoch */
static int xepoch; /* last second epoch */
static int zepoch; /* last run epoch */
static int zcount; /* last run end time */
static int scount; /* seconds counter */
static int syncnt; /* run length counter */
static int maxrun; /* longest run length */
static int mepoch; /* longest run end epoch */
static int mcount; /* longest run end time */
static int avgcnt; /* averaging interval counter */
static int avginc; /* averaging ratchet */
static int iniflg; /* initialization flag */
char tbuf[TBUF]; /* monitor buffer */
double dtemp;
int tmp2;
pp = peer->procptr;
up = pp->unitptr;
if (!iniflg) {
iniflg = 1;
ZERO(epoch_mf);
}
/*
* If the signal amplitude or SNR fall below thresholds, dim the
* second sync lamp and wait for hotter ions. If no stations are
* heard, we are either in a probe cycle or the ions are really
* cold.
*/
scount++;
if (up->epomax < STHR || up->eposnr < SSNR) {
up->status &= ~(SSYNC | FGATE);
avgcnt = syncnt = maxrun = 0;
return;
}
if (!(up->status & (SELV | SELH)))
return;
/*
* A three-stage median filter is used to help denoise the
* second sync pulse. The median sample becomes the candidate
* epoch.
*/
epoch_mf[2] = epoch_mf[1];
epoch_mf[1] = epoch_mf[0];
epoch_mf[0] = epopos;
if (epoch_mf[0] > epoch_mf[1]) {
if (epoch_mf[1] > epoch_mf[2])
tepoch = epoch_mf[1]; /* 0 1 2 */
else if (epoch_mf[2] > epoch_mf[0])
tepoch = epoch_mf[0]; /* 2 0 1 */
else
tepoch = epoch_mf[2]; /* 0 2 1 */
} else {
if (epoch_mf[1] < epoch_mf[2])
tepoch = epoch_mf[1]; /* 2 1 0 */
else if (epoch_mf[2] < epoch_mf[0])
tepoch = epoch_mf[0]; /* 1 0 2 */
else
tepoch = epoch_mf[2]; /* 1 2 0 */
}
/*
* If the epoch candidate is the same as the last one, increment
* the run counter. If not, save the length, epoch and end
* time of the current run for use later and reset the counter.
* The epoch is considered valid if the run is at least SCMP
* (10) s, the minute is synchronized and the interval since the
* last epoch is not greater than the averaging interval. Thus,
* after a long absence, the program will wait a full averaging
* interval while the comb filter charges up and noise
* dissapates..
*/
tmp2 = (tepoch - xepoch) % WWV_SEC;
if (tmp2 == 0) {
syncnt++;
if (syncnt > SCMP && up->status & MSYNC && (up->status &
FGATE || scount - zcount <= up->avgint)) {
up->status |= SSYNC;
up->yepoch = tepoch;
}
} else if (syncnt >= maxrun) {
maxrun = syncnt;
mcount = scount;
mepoch = xepoch;
syncnt = 0;
}
if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status &
MSYNC)) {
snprintf(tbuf, sizeof(tbuf),
"wwv1 %04x %3d %4d %5.0f %5.1f %5d %4d %4d %4d",
up->status, up->gain, tepoch, up->epomax,
up->eposnr, tmp2, avgcnt, syncnt,
maxrun);
record_clock_stats(&peer->srcadr, tbuf);
#ifdef DEBUG
if (debug)
printf("%s\n", tbuf);
#endif /* DEBUG */
}
avgcnt++;
if (avgcnt < up->avgint) {
xepoch = tepoch;
return;
}
/*
* The sample clock frequency is disciplined using a first-order
* feedback loop with time constant consistent with the Allan
* intercept of typical computer clocks. During each averaging
* interval the candidate epoch at the end of the longest run is
* determined. If the longest run is zero, all epoches in the
* interval are different, so the candidate epoch is the current
* epoch. The frequency update is computed from the candidate
* epoch difference (125-us units) and time difference (seconds)
* between updates.
*/
if (syncnt >= maxrun) {
maxrun = syncnt;
mcount = scount;
mepoch = xepoch;
}
xepoch = tepoch;
if (maxrun == 0) {
mepoch = tepoch;
mcount = scount;
}
/*
* The master clock runs at the codec sample frequency of 8000
* Hz, so the intrinsic time resolution is 125 us. The frequency
* resolution ranges from 18 PPM at the minimum averaging
* interval of 8 s to 0.12 PPM at the maximum interval of 1024
* s. An offset update is determined at the end of the longest
* run in each averaging interval. The frequency adjustment is
* computed from the difference between offset updates and the
* interval between them.
*
* The maximum frequency adjustment ranges from 187 PPM at the
* minimum interval to 1.5 PPM at the maximum. If the adjustment
* exceeds the maximum, the update is discarded and the
* hysteresis counter is decremented. Otherwise, the frequency
* is incremented by the adjustment, but clamped to the maximum
* 187.5 PPM. If the update is less than half the maximum, the
* hysteresis counter is incremented. If the counter increments
* to +3, the averaging interval is doubled and the counter set
* to zero; if it decrements to -3, the interval is halved and
* the counter set to zero.
*/
dtemp = (mepoch - zepoch) % WWV_SEC;
if (up->status & FGATE) {
if (fabs(dtemp) < MAXFREQ * MINAVG) {
up->freq += (dtemp / 2.) / ((mcount - zcount) *
FCONST);
if (up->freq > MAXFREQ)
up->freq = MAXFREQ;
else if (up->freq < -MAXFREQ)
up->freq = -MAXFREQ;
if (fabs(dtemp) < MAXFREQ * MINAVG / 2.) {
if (avginc < 3) {
avginc++;
} else {
if (up->avgint < MAXAVG) {
up->avgint <<= 1;
avginc = 0;
}
}
}
} else {
if (avginc > -3) {
avginc--;
} else {
if (up->avgint > MINAVG) {
up->avgint >>= 1;
avginc = 0;
}
}
}
}
if (pp->sloppyclockflag & CLK_FLAG4) {
snprintf(tbuf, sizeof(tbuf),
"wwv2 %04x %5.0f %5.1f %5d %4d %4d %4d %4.0f %7.2f",
up->status, up->epomax, up->eposnr, mepoch,
up->avgint, maxrun, mcount - zcount, dtemp,
up->freq * 1e6 / WWV_SEC);
record_clock_stats(&peer->srcadr, tbuf);
#ifdef DEBUG
if (debug)
printf("%s\n", tbuf);
#endif /* DEBUG */
}
/*
* This is a valid update; set up for the next interval.
*/
up->status |= FGATE;
zepoch = mepoch;
zcount = mcount;
avgcnt = syncnt = maxrun = 0;
}
/*
* wwv_epoch - epoch scanner
*
* This routine extracts data signals from the 100-Hz subcarrier. It
* scans the receiver second epoch to determine the signal amplitudes
* and pulse timings. Receiver synchronization is determined by the
* minute sync pulse detected in the wwv_rf() routine and the second
* sync pulse detected in the wwv_epoch() routine. The transmitted
* signals are delayed by the propagation delay, receiver delay and
* filter delay of this program. Delay corrections are introduced
* separately for WWV and WWVH.
*
* Most communications radios use a highpass filter in the audio stages,
* which can do nasty things to the subcarrier phase relative to the
* sync pulses. Therefore, the data subcarrier reference phase is
* disciplined using the hardlimited quadrature-phase signal sampled at
* the same time as the in-phase signal. The phase tracking loop uses
* phase adjustments of plus-minus one sample (125 us).
*/
static void
wwv_epoch(
struct peer *peer /* peer structure pointer */
)
{
struct refclockproc *pp;
struct wwvunit *up;
struct chan *cp;
static double sigmin, sigzer, sigone, engmax, engmin;
pp = peer->procptr;
up = pp->unitptr;
/*
* Find the maximum minute sync pulse energy for both the
* WWV and WWVH stations. This will be used later for channel
* and station mitigation. Also set the seconds epoch at 800 ms
* well before the end of the second to make sure we never set
* the epoch backwards.
*/
cp = &up->mitig[up->achan];
if (cp->wwv.amp > cp->wwv.syneng)
cp->wwv.syneng = cp->wwv.amp;
if (cp->wwvh.amp > cp->wwvh.syneng)
cp->wwvh.syneng = cp->wwvh.amp;
if (up->rphase == 800 * MS)
up->repoch = up->yepoch;
/*
* Use the signal amplitude at epoch 15 ms as the noise floor.
* This gives a guard time of +-15 ms from the beginning of the
* second until the second pulse rises at 30 ms. There is a
* compromise here; we want to delay the sample as long as
* possible to give the radio time to change frequency and the
* AGC to stabilize, but as early as possible if the second
* epoch is not exact.
*/
if (up->rphase == 15 * MS)
sigmin = sigzer = sigone = up->irig;
/*
* Latch the data signal at 200 ms. Keep this around until the
* end of the second. Use the signal energy as the peak to
* compute the SNR. Use the Q sample to adjust the 100-Hz
* reference oscillator phase.
*/
if (up->rphase == 200 * MS) {
sigzer = up->irig;
engmax = sqrt(up->irig * up->irig + up->qrig *
up->qrig);
up->datpha = up->qrig / up->avgint;
if (up->datpha >= 0) {
up->datapt++;
if (up->datapt >= 80)
up->datapt -= 80;
} else {
up->datapt--;
if (up->datapt < 0)
up->datapt += 80;
}
}
/*
* Latch the data signal at 500 ms. Keep this around until the
* end of the second.
*/
else if (up->rphase == 500 * MS)
sigone = up->irig;
/*
* At the end of the second crank the clock state machine and
* adjust the codec gain. Note the epoch is buffered from the
* center of the second in order to avoid jitter while the
* seconds synch is diddling the epoch. Then, determine the true
* offset and update the median filter in the driver interface.
*
* Use the energy at the end of the second as the noise to
* compute the SNR for the data pulse. This gives a better
* measurement than the beginning of the second, especially when
* returning from the probe channel. This gives a guard time of
* 30 ms from the decay of the longest pulse to the rise of the
* next pulse.
*/
up->rphase++;
if (up->mphase % WWV_SEC == up->repoch) {
up->status &= ~(DGATE | BGATE);
engmin = sqrt(up->irig * up->irig + up->qrig *
up->qrig);
up->datsig = engmax;
up->datsnr = wwv_snr(engmax, engmin);
/*
* If the amplitude or SNR is below threshold, average a
* 0 in the the integrators; otherwise, average the
* bipolar signal. This is done to avoid noise polution.
*/
if (engmax < DTHR || up->datsnr < DSNR) {
up->status |= DGATE;
wwv_rsec(peer, 0);
} else {
sigzer -= sigone;
sigone -= sigmin;
wwv_rsec(peer, sigone - sigzer);
}
if (up->status & (DGATE | BGATE))
up->errcnt++;
if (up->errcnt > MAXERR)
up->alarm |= LOWERR;
wwv_gain(peer);
cp = &up->mitig[up->achan];
cp->wwv.syneng = 0;
cp->wwvh.syneng = 0;
up->rphase = 0;
}
}
/*
* wwv_rsec - process receiver second
*
* This routine is called at the end of each receiver second to
* implement the per-second state machine. The machine assembles BCD
* digit bits, decodes miscellaneous bits and dances the leap seconds.
*
* Normally, the minute has 60 seconds numbered 0-59. If the leap
* warning bit is set, the last minute (1439) of 30 June (day 181 or 182
* for leap years) or 31 December (day 365 or 366 for leap years) is
* augmented by one second numbered 60. This is accomplished by
* extending the minute interval by one second and teaching the state
* machine to ignore it.
*/
static void
wwv_rsec(
struct peer *peer, /* peer structure pointer */
double bit
)
{
static int iniflg; /* initialization flag */
static double bcddld[4]; /* BCD data bits */
static double bitvec[61]; /* bit integrator for misc bits */
struct refclockproc *pp;
struct wwvunit *up;
struct chan *cp;
struct sync *sp, *rp;
char tbuf[TBUF]; /* monitor buffer */
int sw, arg, nsec;
pp = peer->procptr;
up = pp->unitptr;
if (!iniflg) {
iniflg = 1;
ZERO(bitvec);
}
/*
* The bit represents the probability of a hit on zero (negative
* values), a hit on one (positive values) or a miss (zero
* value). The likelihood vector is the exponential average of
* these probabilities. Only the bits of this vector
* corresponding to the miscellaneous bits of the timecode are
* used, but it's easier to do them all. After that, crank the
* seconds state machine.
*/
nsec = up->rsec;
up->rsec++;
bitvec[nsec] += (bit - bitvec[nsec]) / TCONST;
sw = progx[nsec].sw;
arg = progx[nsec].arg;
/*
* The minute state machine. Fly off to a particular section as
* directed by the transition matrix and second number.
*/
switch (sw) {
/*
* Ignore this second.
*/
case IDLE: /* 9, 45-49 */
break;
/*
* Probe channel stuff
*
* The WWV/H format contains data pulses in second 59 (position
* identifier) and second 1, but not in second 0. The minute
* sync pulse is contained in second 0. At the end of second 58
* QSY to the probe channel, which rotates in turn over all
* WWV/H frequencies. At the end of second 0 measure the minute
* sync pulse. At the end of second 1 measure the data pulse and
* QSY back to the data channel. Note that the actions commented
* here happen at the end of the second numbered as shown.
*
* At the end of second 0 save the minute sync amplitude latched
* at 800 ms as the signal later used to calculate the SNR.
*/
case SYNC2: /* 0 */
cp = &up->mitig[up->achan];
cp->wwv.synmax = cp->wwv.syneng;
cp->wwvh.synmax = cp->wwvh.syneng;
break;
/*
* At the end of second 1 use the minute sync amplitude latched
* at 800 ms as the noise to calculate the SNR. If the minute
* sync pulse and SNR are above thresholds and the data pulse
* amplitude and SNR are above thresolds, shift a 1 into the
* station reachability register; otherwise, shift a 0. The
* number of 1 bits in the last six intervals is a component of
* the channel metric computed by the wwv_metric() routine.
* Finally, QSY back to the data channel.
*/
case SYNC3: /* 1 */
cp = &up->mitig[up->achan];
/*
* If synchronized to a station, restart if no stations
* have been heard within the PANIC timeout (2 days). If
* not and the minute digit has been found, restart if
* not synchronized withing the SYNCH timeout (40 m). If
* not, restart if the unit digit has not been found
* within the DATA timeout (15 m).
*/
if (up->status & INSYNC) {
if (up->watch > PANIC) {
wwv_newgame(peer);
return;
}
} else if (up->status & DSYNC) {
if (up->watch > SYNCH) {
wwv_newgame(peer);
return;
}
} else if (up->watch > DATA) {
wwv_newgame(peer);
return;
}
wwv_newchan(peer);
break;
/*
* Save the bit probability in the BCD data vector at the index
* given by the argument. Bits not used in the digit are forced
* to zero.
*/
case COEF1: /* 4-7 */
bcddld[arg] = bit;
break;
/*
* Correlate coefficient vector with each valid digit vector and
* save in decoding matrix. We step through the decoding matrix
* digits correlating each with the coefficients and saving the
* greatest and the next lower for later SNR calculation.
*/
case DECIM2: /* 29 */
wwv_corr4(peer, &up->decvec[arg], bcddld, bcd2);
break;
case DECIM3: /* 44 */
wwv_corr4(peer, &up->decvec[arg], bcddld, bcd3);
break;
case DECIM6: /* 19 */
wwv_corr4(peer, &up->decvec[arg], bcddld, bcd6);
break;
/*
* Miscellaneous bits. If above the positive threshold, declare
* 1; if below the negative threshold, declare 0; otherwise
* raise the BGATE bit. The design is intended to avoid
* integrating noise under low SNR conditions.
*/
case MSC20: /* 55 */
wwv_corr4(peer, &up->decvec[YR + 1], bcddld, bcd9);
/* fall through */
case MSCBIT: /* 2-3, 50, 56-57 */
if (bitvec[nsec] > BTHR) {
if (!(up->misc & arg))
up->alarm |= CMPERR;
up->misc |= arg;
} else if (bitvec[nsec] < -BTHR) {
if (up->misc & arg)
up->alarm |= CMPERR;
up->misc &= ~arg;
} else {
up->status |= BGATE;
}
break;
/*
* Save the data channel gain, then QSY to the probe channel and
* dim the seconds comb filters. The www_newchan() routine will
* light them back up.
*/
case MSC21: /* 58 */
if (bitvec[nsec] > BTHR) {
if (!(up->misc & arg))
up->alarm |= CMPERR;
up->misc |= arg;
} else if (bitvec[nsec] < -BTHR) {
if (up->misc & arg)
up->alarm |= CMPERR;
up->misc &= ~arg;
} else {
up->status |= BGATE;
}
up->status &= ~(SELV | SELH);
#ifdef ICOM
if (up->fd_icom > 0) {
up->schan = (up->schan + 1) % NCHAN;
wwv_qsy(peer, up->schan);
} else {
up->mitig[up->achan].gain = up->gain;
}
#else
up->mitig[up->achan].gain = up->gain;
#endif /* ICOM */
break;
/*
* The endgames
*
* During second 59 the receiver and codec AGC are settling
* down, so the data pulse is unusable as quality metric. If
* LEPSEC is set on the last minute of 30 June or 31 December,
* the transmitter and receiver insert an extra second (60) in
* the timescale and the minute sync repeats the second. Once
* leaps occurred at intervals of about 18 months, but the last
* leap before the most recent leap in 1995 was in 1998.
*/
case MIN1: /* 59 */
if (up->status & LEPSEC)
break;
/*
* The radio clock is set if the alarm bits are all zero. After that,
* the time is considered valid if the second sync bit is lit. It should
* not be a surprise, especially if the radio is not tunable, that
* sometimes no stations are above the noise and the integrators
* discharge below the thresholds. We assume that, after a day of signal
* loss, the minute sync epoch will be in the same second. This requires
* the codec frequency be accurate within 6 PPM. Practical experience
* shows the frequency typically within 0.1 PPM, so after a day of
* signal loss, the time should be within 8.6 ms..
*/
static void
wwv_clock(
struct peer *peer /* peer unit pointer */
)
{
struct refclockproc *pp;
struct wwvunit *up;
l_fp offset; /* offset in NTP seconds */
/*
* wwv_corr4 - determine maximum-likelihood digit
*
* This routine correlates the received digit vector with the BCD
* coefficient vectors corresponding to all valid digits at the given
* position in the decoding matrix. The maximum value corresponds to the
* maximum-likelihood digit, while the ratio of this value to the next
* lower value determines the likelihood function. Note that, if the
* digit is invalid, the likelihood vector is averaged toward a miss.
*/
static void
wwv_corr4(
struct peer *peer, /* peer unit pointer */
struct decvec *vp, /* decoding table pointer */
double data[], /* received data vector */
double tab[][4] /* correlation vector array */
)
{
struct refclockproc *pp;
struct wwvunit *up;
double topmax, nxtmax; /* metrics */
double acc; /* accumulator */
char tbuf[TBUF]; /* monitor buffer */
int mldigit; /* max likelihood digit */
int i, j;
pp = peer->procptr;
up = pp->unitptr;
/*
* Correlate digit vector with each BCD coefficient vector. If
* any BCD digit bit is bad, consider all bits a miss. Until the
* minute units digit has been resolved, don't to anything else.
* Note the SNR is calculated as the ratio of the largest
* likelihood value to the next largest likelihood value.
*/
mldigit = 0;
topmax = nxtmax = -MAXAMP;
for (i = 0; tab[i][0] != 0; i++) {
acc = 0;
for (j = 0; j < 4; j++)
acc += data[j] * tab[i][j];
acc = (vp->like[i] += (acc - vp->like[i]) / TCONST);
if (acc > topmax) {
nxtmax = topmax;
topmax = acc;
mldigit = i;
} else if (acc > nxtmax) {
nxtmax = acc;
}
}
vp->digprb = topmax;
vp->digsnr = wwv_snr(topmax, nxtmax);
/*
* The current maximum-likelihood digit is compared to the last
* maximum-likelihood digit. If different, the compare counter
* and maximum-likelihood digit are reset. When the compare
* counter reaches the BCMP threshold (3), the digit is assumed
* correct. When the compare counter of all nine digits have
* reached threshold, the clock is assumed correct.
*
* Note that the clock display digit is set before the compare
* counter has reached threshold; however, the clock display is
* not considered correct until all nine clock digits have
* reached threshold. This is intended as eye candy, but avoids
* mistakes when the signal is low and the SNR is very marginal.
*/
if (vp->digprb < BTHR || vp->digsnr < BSNR) {
up->status |= BGATE;
} else {
if (vp->digit != mldigit) {
up->alarm |= CMPERR;
if (vp->count > 0)
vp->count--;
if (vp->count == 0)
vp->digit = mldigit;
} else {
if (vp->count < BCMP)
vp->count++;
if (vp->count == BCMP) {
up->status |= DSYNC;
up->digcnt++;
}
}
}
if ((pp->sloppyclockflag & CLK_FLAG4) && !(up->status &
INSYNC)) {
snprintf(tbuf, sizeof(tbuf),
"wwv4 %2d %04x %3d %4d %5.0f %2d %d %d %d %5.0f %5.1f",
up->rsec - 1, up->status, up->gain, up->yepoch,
up->epomax, vp->radix, vp->digit, mldigit,
vp->count, vp->digprb, vp->digsnr);
record_clock_stats(&peer->srcadr, tbuf);
#ifdef DEBUG
if (debug)
printf("%s\n", tbuf);
#endif /* DEBUG */
}
}
/*
* wwv_tsec - transmitter minute processing
*
* This routine is called at the end of the transmitter minute. It
* implements a state machine that advances the logical clock subject to
* the funny rules that govern the conventional clock and calendar.
*/
static void
wwv_tsec(
struct peer *peer /* driver structure pointer */
)
{
struct refclockproc *pp;
struct wwvunit *up;
int minute, day, isleap;
int temp;
pp = peer->procptr;
up = pp->unitptr;
/*
* Advance minute unit of the day. Don't propagate carries until
* the unit minute digit has been found.
*/
temp = carry(&up->decvec[MN]); /* minute units */
if (!(up->status & DSYNC))
return;
/*
* Propagate carries through the day.
*/
if (temp == 0) /* carry minutes */
temp = carry(&up->decvec[MN + 1]);
if (temp == 0) /* carry hours */
temp = carry(&up->decvec[HR]);
if (temp == 0)
temp = carry(&up->decvec[HR + 1]);
// XXX: Does temp have an expected value here?
/*
* Decode the current minute and day. Set leap day if the
* timecode leap bit is set on 30 June or 31 December. Set leap
* minute if the last minute on leap day, but only if the clock
* is syncrhronized. This code fails in 2400 AD.
*/
minute = up->decvec[MN].digit + up->decvec[MN + 1].digit *
10 + up->decvec[HR].digit * 60 + up->decvec[HR +
1].digit * 600;
day = up->decvec[DA].digit + up->decvec[DA + 1].digit * 10 +
up->decvec[DA + 2].digit * 100;
/*
* Set the leap bit on the last minute of the leap day.
*/
isleap = up->decvec[YR].digit & 0x3;
if (up->misc & SECWAR && up->status & INSYNC) {
if ((day == (isleap ? 182 : 183) || day == (isleap ?
365 : 366)) && minute == 1439)
up->status |= LEPSEC;
}
/*
* Roll the day if this the first minute and propagate carries
* through the year.
*/
if (minute != 1440)
return;
// minute = 0;
while (carry(&up->decvec[HR]) != 0); /* advance to minute 0 */
while (carry(&up->decvec[HR + 1]) != 0);
day++;
temp = carry(&up->decvec[DA]); /* carry days */
if (temp == 0)
temp = carry(&up->decvec[DA + 1]);
if (temp == 0)
temp = carry(&up->decvec[DA + 2]);
// XXX: Is there an expected value of temp here?
/*
* Roll the year if this the first day and propagate carries
* through the century.
*/
if (day != (isleap ? 365 : 366))
return;
// day = 1;
while (carry(&up->decvec[DA]) != 1); /* advance to day 1 */
while (carry(&up->decvec[DA + 1]) != 0);
while (carry(&up->decvec[DA + 2]) != 0);
temp = carry(&up->decvec[YR]); /* carry years */
if (temp == 0)
carry(&up->decvec[YR + 1]);
}
/*
* carry - process digit
*
* This routine rotates a likelihood vector one position and increments
* the clock digit modulo the radix. It returns the new clock digit or
* zero if a carry occurred. Once synchronized, the clock digit will
* match the maximum-likelihood digit corresponding to that position.
*/
static int
carry(
struct decvec *dp /* decoding table pointer */
)
{
int temp;
int j;
/*
* wwv_snr - compute SNR or likelihood function
*/
static double
wwv_snr(
double signal, /* signal */
double noise /* noise */
)
{
double rval;
/*
* This is a little tricky. Due to the way things are measured,
* either or both the signal or noise amplitude can be negative
* or zero. The intent is that, if the signal is negative or
* zero, the SNR must always be zero. This can happen with the
* subcarrier SNR before the phase has been aligned. On the
* other hand, in the likelihood function the "noise" is the
* next maximum down from the peak and this could be negative.
* However, in this case the SNR is truly stupendous, so we
* simply cap at MAXSNR dB (40).
*/
if (signal <= 0) {
rval = 0;
} else if (noise <= 0) {
rval = MAXSNR;
} else {
rval = 20. * log10(signal / noise);
if (rval > MAXSNR)
rval = MAXSNR;
}
return (rval);
}
/*
* wwv_newchan - change to new data channel
*
* The radio actually appears to have ten channels, one channel for each
* of five frequencies and each of two stations (WWV and WWVH), although
* if not tunable only the DCHAN channel appears live. While the radio
* is tuned to the working data channel frequency and station for most
* of the minute, during seconds 59, 0 and 1 the radio is tuned to a
* probe frequency in order to search for minute sync pulse and data
* subcarrier from other transmitters.
*
* The search for WWV and WWVH operates simultaneously, with WWV minute
* sync pulse at 1000 Hz and WWVH at 1200 Hz. The probe frequency
* rotates each minute over 2.5, 5, 10, 15 and 20 MHz in order and yes,
* we all know WWVH is dark on 20 MHz, but few remember when WWV was lit
* on 25 MHz.
*
* This routine selects the best channel using a metric computed from
* the reachability register and minute pulse amplitude. Normally, the
* award goes to the the channel with the highest metric; but, in case
* of ties, the award goes to the channel with the highest minute sync
* pulse amplitude and then to the highest frequency.
*
* The routine performs an important squelch function to keep dirty data
* from polluting the integrators. In order to consider a station valid,
* the metric must be at least MTHR (13); otherwise, the station select
* bits are cleared so the second sync is disabled and the data bit
* integrators averaged to a miss.
*/
static int
wwv_newchan(
struct peer *peer /* peer structure pointer */
)
{
struct refclockproc *pp;
struct wwvunit *up;
struct sync *sp, *rp;
double rank, dtemp;
int i, j, rval;
pp = peer->procptr;
up = pp->unitptr;
/*
* Search all five station pairs looking for the channel with
* maximum metric.
*/
sp = NULL;
j = 0;
rank = 0;
for (i = 0; i < NCHAN; i++) {
rp = &up->mitig[i].wwvh;
dtemp = rp->metric;
if (dtemp >= rank) {
rank = dtemp;
sp = rp;
j = i;
}
rp = &up->mitig[i].wwv;
dtemp = rp->metric;
if (dtemp >= rank) {
rank = dtemp;
sp = rp;
j = i;
}
}
/*
* If the strongest signal is less than the MTHR threshold (13),
* we are beneath the waves, so squelch the second sync and
* advance to the next station. This makes sure all stations are
* scanned when the ions grow dim. If the strongest signal is
* greater than the threshold, tune to that frequency and
* transmitter QTH.
*/
up->status &= ~(SELV | SELH);
if (rank < MTHR) {
up->dchan = (up->dchan + 1) % NCHAN;
if (up->status & METRIC) {
up->status &= ~METRIC;
refclock_report(peer, CEVNT_PROP);
}
rval = FALSE;
} else {
up->dchan = j;
up->sptr = sp;
memcpy(&pp->refid, sp->refid, 4);
peer->refid = pp->refid;
up->status |= METRIC;
if (sp->select & SELV) {
up->status |= SELV;
up->pdelay = pp->fudgetime1;
} else if (sp->select & SELH) {
up->status |= SELH;
up->pdelay = pp->fudgetime2;
} else {
up->pdelay = 0;
}
rval = TRUE;
}
#ifdef ICOM
if (up->fd_icom > 0)
wwv_qsy(peer, up->dchan);
#endif /* ICOM */
return (rval);
}
/*
* wwv_newgame - reset and start over
*
* There are three conditions resulting in a new game:
*
* 1 After finding the minute pulse (MSYNC lit), going 15 minutes
* (DATA) without finding the unit seconds digit.
*
* 2 After finding good data (DSYNC lit), going more than 40 minutes
* (SYNCH) without finding station sync (INSYNC lit).
*
* 3 After finding station sync (INSYNC lit), going more than 2 days
* (PANIC) without finding any station.
*/
static void
wwv_newgame(
struct peer *peer /* peer structure pointer */
)
{
struct refclockproc *pp;
struct wwvunit *up;
struct chan *cp;
int i;
pp = peer->procptr;
up = pp->unitptr;
/*
* Initialize strategic values. Note we set the leap bits
* NOTINSYNC and the refid "NONE".
*/
if (up->status)
up->errflg = CEVNT_TIMEOUT;
peer->leap = LEAP_NOTINSYNC;
up->watch = up->status = up->alarm = 0;
up->avgint = MINAVG;
up->freq = 0;
up->gain = MAXGAIN / 2;
/*
* Initialize the station processes for audio gain, select bit,
* station/frequency identifier and reference identifier. Start
* probing at the strongest channel or the default channel if
* nothing heard.
*/
memset(up->mitig, 0, sizeof(up->mitig));
for (i = 0; i < NCHAN; i++) {
cp = &up->mitig[i];
cp->gain = up->gain;
cp->wwv.select = SELV;
snprintf(cp->wwv.refid, sizeof(cp->wwv.refid), "WV%.0f",
floor(qsy[i]));
cp->wwvh.select = SELH;
snprintf(cp->wwvh.refid, sizeof(cp->wwvh.refid), "WH%.0f",
floor(qsy[i]));
}
up->dchan = (DCHAN + NCHAN - 1) % NCHAN;
wwv_newchan(peer);
up->schan = up->dchan;
}
/*
* wwv_metric - compute station metric
*
* The most significant bits represent the number of ones in the
* station reachability register. The least significant bits represent
* the minute sync pulse amplitude. The combined value is scaled 0-100.
*/
double
wwv_metric(
struct sync *sp /* station pointer */
)
{
double dtemp;
#ifdef ICOM
/*
* wwv_qsy - Tune ICOM receiver
*
* This routine saves the AGC for the current channel, switches to a new
* channel and restores the AGC for that channel. If a tunable receiver
* is not available, just fake it.
*/
static int
wwv_qsy(
struct peer *peer, /* peer structure pointer */
int chan /* channel */
)
{
int rval = 0;
struct refclockproc *pp;
struct wwvunit *up;
/*
* timecode - assemble timecode string and length
*
* Prettytime format - similar to Spectracom
*
* sq yy ddd hh:mm:ss ld dut lset agc iden sig errs freq avgt
*
* s sync indicator ('?' or ' ')
* q error bits (hex 0-F)
* yyyy year of century
* ddd day of year
* hh hour of day
* mm minute of hour
* ss second of minute)
* l leap second warning (' ' or 'L')
* d DST state ('S', 'D', 'I', or 'O')
* dut DUT sign and magnitude (0.1 s)
* lset minutes since last clock update
* agc audio gain (0-255)
* iden reference identifier (station and frequency)
* sig signal quality (0-100)
* errs bit errors in last minute
* freq frequency offset (PPM)
* avgt averaging time (s)
*/
static int
timecode(
struct wwvunit *up, /* driver structure pointer */
char * tc, /* target string */
size_t tcsiz /* target max chars */
)
{
struct sync *sp;
int year, day, hour, minute, second, dut;
char synchar, leapchar, dst;
char cptr[50];
/*
* wwv_gain - adjust codec gain
*
* This routine is called at the end of each second. During the second
* the number of signal clips above the MAXAMP threshold (6000). If
* there are no clips, the gain is bumped up; if there are more than
* MAXCLP clips (100), it is bumped down. The decoder is relatively
* insensitive to amplitude, so this crudity works just peachy. The
* routine also jiggles the input port and selectively mutes the
* monitor.
*/
static void
wwv_gain(
struct peer *peer /* peer structure pointer */
)
{
struct refclockproc *pp;
struct wwvunit *up;
pp = peer->procptr;
up = pp->unitptr;
/*
* Apparently, the codec uses only the high order bits of the
* gain control field. Thus, it may take awhile for changes to
* wiggle the hardware bits.
*/
if (up->clipcnt == 0) {
up->gain += 4;
if (up->gain > MAXGAIN)
up->gain = MAXGAIN;
} else if (up->clipcnt > MAXCLP) {
up->gain -= 4;
if (up->gain < 0)
up->gain = 0;
}
audio_gain(up->gain, up->mongain, up->port);
up->clipcnt = 0;
#if DEBUG
if (debug > 1)
audio_show();
#endif
}
