/* * refclock_wwv - clock driver for NIST WWV/H time/frequency station */ #ifdef HAVE_CONFIG_H #include <config.h> #endif #if defined(REFCLOCK) && defined(CLOCK_WWV) #include "ntpd.h" #include "ntp_io.h" #include "ntp_refclock.h" #include "ntp_calendar.h" #include "ntp_stdlib.h" #include "audio.h" #include <stdio.h> #include <ctype.h> #include <math.h> #ifdef HAVE_SYS_IOCTL_H # include <sys/ioctl.h> #endif /* HAVE_SYS_IOCTL_H */ #define ICOM 1 /* undefine to suppress ICOM code */ #ifdef ICOM #include "icom.h" #endif /* ICOM */ /* * 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. Transmikssions are made continuously on * 2.5, 5, 10, 15 and 20 MHz in AM mode. An ordinary 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 day and * night. * * The driver receives, demodulates and decodes the radio signals when * connected to the audio codec of a Sun workstation running SunOS or * Solaris, and with a little help, other workstations with similar * codecs or sound cards. 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.htm. 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. */ /* * Interface definitions */ #define DEVICE_AUDIO "/dev/audio" /* audio device name */ #define PRECISION (-10) /* precision assumed (about 1 ms) */ #define REFID "NONE" /* reference ID */ #define DESCRIPTION "WWV/H Audio Demodulator/Decoder" /* WRU */ #define SECOND 8000 /* second epoch (sample rate) (Hz) */ #define MINUTE (SECOND * 60) /* minute epoch */ #define OFFSET 128 /* companded sample offset */ #define SIZE 256 /* decompanding table size */ #define MAXSIG 6000. /* maximum signal level reference */ #define MAXSNR 30. /* max SNR reference */ #define DGAIN 20. /* data channel gain reference */ #define SGAIN 10. /* sync channel gain reference */ #define MAXFREQ (125e-6 * SECOND) /* freq tolerance (.0125%) */ #define PI 3.1415926535 /* the real thing */ #define DATSIZ (170 * MS) /* data matched filter size */ #define SYNSIZ (800 * MS) /* minute sync matched filter size */ #define UTCYEAR 72 /* the first UTC year */ #define MAXERR 30 /* max data bit errors in minute */ #define NCHAN 5 /* number of channels */ /* * Macroni */ #define MOD(x, y) ((x) < 0 ? -(-(x) % (y)) : (x) % (y)) /* * General purpose status bits (status) * * Notes: SELV and/or SELH are set when the minute sync pulse from * either or both WWV and/or WWVH stations has been heard. MSYNC is set * when the minute sync pulse has been acquired and never reset. SSYNC * is set when the second sync pulse has been acquired and cleared by * watchdog or signal loss. DSYNC is set when the minutes unit digit has * reached the threshold and INSYNC is set when if all nine digits have * reached the threshold and never cleared. * * DGATE is set if a data bit is invalid, BGATE is set if a BCD digit * bit is invalid. SFLAG is set when during seconds 59, 0 and 1 while * probing for alternate frequencies. LEPSEC is set when the SECWAR of * the timecode is set on the last second of 30 June or 31 December. At * the end of this minute both the receiver and transmitter insert * second 60 in the minute 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 DGATE 0x0010 /* data bit error */ #define BGATE 0x0020 /* BCD digit bit error */ #define SFLAG 0x1000 /* probe flag */ #define LEPSEC 0x2000 /* leap second in progress */ /* * Station scoreboard bits (select) * * These are used to establish the signal quality for each of the five * frequencies and two stations. */ #define JITRNG 0x0001 /* jitter above threshold */ #define SYNCNG 0x0002 /* sync below threshold or SNR */ #define DATANG 0x0004 /* data below threshold or SNR */ #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. There are four * four-bit nibble fields in the word, each corresponding to a specific * alarm condition. At the end of each second, the word is shifted left * one position and the least significant bit of each nibble cleared. * This bit can be set during the next minute if the associated alarm * condition is raised. This provides a way to remember alarm conditions * up to four minutes. * * If not tracking both minute sync and second sync, the SYNERR alarm is * raised. The data error counter is incremented for each invalid data * bit. If too many data bit errors are encountered in one minute, the * MODERR alarm is raised. The DECERR alarm is raised if a maximum * likelihood digit fails to compare with the current clock digit. If * the probability of any miscellaneous bit or any digit falls below the * threshold, the SYMERR alarm is raised. */ #define DECERR 0 /* BCD digit compare error */ #define SYMERR 4 /* low bit or digit probability */ #define MODERR 8 /* too many data bit errors */ #define SYNERR 12 /* not synchronized to station */ /* * Watchdog 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 5 /* acquisition timeout */ #define HSPEC 15 /* second sync timeout */ #define DIGIT 30 /* minute unit digit timeout */ #define PANIC (4 * 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 receiver. The values defined here may be on the * adventurous side in the interest of the highest sensitivity. */ #define ATHR 2000 /* acquisition amplitude threshold */ #define ASNR 6.0 /* acquisition SNR threshold (dB) */ #define AWND 50 /* acquisition window threshold (ms) */ #define AMIN 3 /* acquisition min compare count */ #define AMAX 6 /* max compare count */ #define QTHR 2000 /* QSY amplitude threshold */ #define QSNR 20.0 /* QSY SNR threshold (dB) */ #define STHR 500 /* second sync amplitude threshold */ #define SCMP 10 /* second sync compare threshold */ #define DTHR 1000 /* bit amplitude threshold */ #define DSNR 10.0 /* bit SNR threshold (dB) */ #define BTHR 1000 /* digit probability threshold */ #define BSNR 3.0 /* digit likelihood threshold (dB) */ #define BCMP 5 /* digit compare threshold (dB) */ /* * Tone frequency definitions. */ #define MS 8 /* samples per millisecond */ #define IN100 1 /* 100 Hz 4.5-deg sin table */ #define IN1000 10 /* 1000 Hz 4.5-deg sin table */ #define IN1200 12 /* 1200 Hz 4.5-deg sin table */ /* * Acquisition and tracking time constants. Usually powers of 2. */ #define MINAVG 8 /* min time constant (s) */ #define MAXAVG 7 /* max time constant (log2 s) */ #define TCONST 16 /* minute time constant (s) */ #define SYNCTC (1024 / (1 << MAXAVG)) /* FLL constant (s) */ /* * 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 DST in progress */ #define DST2 0x20 /* 2 DST2 DST change warning */ #define SECWAR 0x40 /* 3 leap second warning */ /* * The total system delay with the DSP93 program is at 22.5 ms, * including the propagation delay from Ft. Collins, CO, to Newark, DE * (8.9 ms), the communications receiver delay and the delay of the * DSP93 program itself. The DSP93 program delay is due mainly to the * 400-Hz FIR bandpass filter (5 ms) and second sync matched filter (5 * ms), leaving about 3.6 ms for the receiver delay and strays. * * The total system delay with this program is estimated at 27.1 ms by * comparison with another PPS-synchronized NTP server over a 10-Mb/s * Ethernet. The propagation and receiver delays are the same as with * the DSP93 program. The program delay is due only to the 600-Hz * IIR bandpass filter (1.1 ms), since other delays have been removed. * Assuming 4.7 ms for the receiver, program and strays, this leaves * 13.5 ms for the audio codec and operating system latencies for a * total of 18.2 ms. as the systematic delay. The additional propagation * delay specific to each receiver location can be programmed in the * fudge time1 and time2 values for WWV and WWVH, respectively. */ #define PDELAY (.0036 + .0011 + .0135) /* net system delay (s) */ /* * Table of sine values at 4.5-degree increments. This is used by the * synchronous matched filter demodulators. The integral of sine-squared * over one complete cycle is PI, so the table is normallized by 1 / PI. */ double sintab[] = { 0.000000e+00, 2.497431e-02, 4.979464e-02, 7.430797e-02, /* 0-3 */ 9.836316e-02, 1.218119e-01, 1.445097e-01, 1.663165e-01, /* 4-7 */ 1.870979e-01, 2.067257e-01, 2.250791e-01, 2.420447e-01, /* 8-11 */ 2.575181e-01, 2.714038e-01, 2.836162e-01, 2.940800e-01, /* 12-15 */ 3.027307e-01, 3.095150e-01, 3.143910e-01, 3.173286e-01, /* 16-19 */ 3.183099e-01, 3.173286e-01, 3.143910e-01, 3.095150e-01, /* 20-23 */ 3.027307e-01, 2.940800e-01, 2.836162e-01, 2.714038e-01, /* 24-27 */ 2.575181e-01, 2.420447e-01, 2.250791e-01, 2.067257e-01, /* 28-31 */ 1.870979e-01, 1.663165e-01, 1.445097e-01, 1.218119e-01, /* 32-35 */ 9.836316e-02, 7.430797e-02, 4.979464e-02, 2.497431e-02, /* 36-39 */ -0.000000e+00, -2.497431e-02, -4.979464e-02, -7.430797e-02, /* 40-43 */ -9.836316e-02, -1.218119e-01, -1.445097e-01, -1.663165e-01, /* 44-47 */ -1.870979e-01, -2.067257e-01, -2.250791e-01, -2.420447e-01, /* 48-51 */ -2.575181e-01, -2.714038e-01, -2.836162e-01, -2.940800e-01, /* 52-55 */ -3.027307e-01, -3.095150e-01, -3.143910e-01, -3.173286e-01, /* 56-59 */ -3.183099e-01, -3.173286e-01, -3.143910e-01, -3.095150e-01, /* 60-63 */ -3.027307e-01, -2.940800e-01, -2.836162e-01, -2.714038e-01, /* 64-67 */ -2.575181e-01, -2.420447e-01, -2.250791e-01, -2.067257e-01, /* 68-71 */ -1.870979e-01, -1.663165e-01, -1.445097e-01, -1.218119e-01, /* 72-75 */ -9.836316e-02, -7.430797e-02, -4.979464e-02, -2.497431e-02, /* 76-79 */ 0.000000e+00}; /* 80 */ /* * 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 conditioned on DSYNC */ #define COEF1 2 /* BCD bit */ #define COEF2 3 /* BCD bit ignored */ #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 /* minute */ #define MIN2 12 /* leap second */ #define SYNC2 13 /* QSY data channel */ #define SYNC3 14 /* QSY data channel */ /* * 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) */ struct progx progx[] = { {SYNC2, 0}, /* 0 latch sync max */ {SYNC3, 0}, /* 1 QSY data channel */ {MSCBIT, DST2}, /* 2 dst2 */ {MSCBIT, SECWAR}, /* 3 lw */ {COEF, 0}, /* 4 1 year units */ {COEF, 1}, /* 5 2 */ {COEF, 2}, /* 6 4 */ {COEF, 3}, /* 7 8 */ {DECIM9, YR}, /* 8 */ {IDLE, 0}, /* 9 p1 */ {COEF1, 0}, /* 10 1 minute units */ {COEF1, 1}, /* 11 2 */ {COEF1, 2}, /* 12 4 */ {COEF1, 3}, /* 13 8 */ {DECIM9, MN}, /* 14 */ {COEF, 0}, /* 15 10 minute tens */ {COEF, 1}, /* 16 20 */ {COEF, 2}, /* 17 40 */ {COEF2, 3}, /* 18 80 (not used) */ {DECIM6, MN + 1}, /* 19 p2 */ {COEF, 0}, /* 20 1 hour units */ {COEF, 1}, /* 21 2 */ {COEF, 2}, /* 22 4 */ {COEF, 3}, /* 23 8 */ {DECIM9, HR}, /* 24 */ {COEF, 0}, /* 25 10 hour tens */ {COEF, 1}, /* 26 20 */ {COEF2, 2}, /* 27 40 (not used) */ {COEF2, 3}, /* 28 80 (not used) */ {DECIM2, HR + 1}, /* 29 p3 */ {COEF, 0}, /* 30 1 day units */ {COEF, 1}, /* 31 2 */ {COEF, 2}, /* 32 4 */ {COEF, 3}, /* 33 8 */ {DECIM9, DA}, /* 34 */ {COEF, 0}, /* 35 10 day tens */ {COEF, 1}, /* 36 20 */ {COEF, 2}, /* 37 40 */ {COEF, 3}, /* 38 80 */ {DECIM9, DA + 1}, /* 39 p4 */ {COEF, 0}, /* 40 100 day hundreds */ {COEF, 1}, /* 41 200 */ {COEF2, 2}, /* 42 400 (not used) */ {COEF2, 3}, /* 43 800 (not used) */ {DECIM3, DA + 2}, /* 44 */ {IDLE, 0}, /* 45 */ {IDLE, 0}, /* 46 */ {IDLE, 0}, /* 47 */ {IDLE, 0}, /* 48 */ {IDLE, 0}, /* 49 p5 */ {MSCBIT, DUTS}, /* 50 dut+- */ {COEF, 0}, /* 51 10 year tens */ {COEF, 1}, /* 52 20 */ {COEF, 2}, /* 53 40 */ {COEF, 3}, /* 54 80 */ {MSC20, DST1}, /* 55 dst1 */ {MSCBIT, DUT1}, /* 56 0.1 dut */ {MSCBIT, DUT2}, /* 57 0.2 */ {MSC21, DUT4}, /* 58 0.4 QSY probe channel */ {MIN1, 0}, /* 59 p6 latch sync min */ {MIN2, 0} /* 60 leap second */ }; /* * BCD coefficients for maximum likelihood digit decode */ #define P15 1. /* max positive number */ #define N15 -1. /* max negative number */ /* * Digits 0-9 */ #define P9 (P15 / 4) /* mark (+1) */ #define N9 (N15 / 4) /* space (-1) */ double bcd9[][4] = { {N9, N9, N9, N9}, /* 0 */ {P9, N9, N9, N9}, /* 1 */ {N9, P9, N9, N9}, /* 2 */ {P9, P9, N9, N9}, /* 3 */ {N9, N9, P9, N9}, /* 4 */ {P9, N9, P9, N9}, /* 5 */ {N9, P9, P9, N9}, /* 6 */ {P9, P9, P9, N9}, /* 7 */ {N9, N9, N9, P9}, /* 8 */ {P9, N9, N9, P9}, /* 9 */ {0, 0, 0, 0} /* backstop */ }; /* * Digits 0-6 (minute tens) */ #define P6 (P15 / 3) /* mark (+1) */ #define N6 (N15 / 3) /* space (-1) */ double bcd6[][4] = { {N6, N6, N6, 0}, /* 0 */ {P6, N6, N6, 0}, /* 1 */ {N6, P6, N6, 0}, /* 2 */ {P6, P6, N6, 0}, /* 3 */ {N6, N6, P6, 0}, /* 4 */ {P6, N6, P6, 0}, /* 5 */ {N6, P6, P6, 0}, /* 6 */ {0, 0, 0, 0} /* backstop */ }; /* * Digits 0-3 (day hundreds) */ #define P3 (P15 / 2) /* mark (+1) */ #define N3 (N15 / 2) /* space (-1) */ double bcd3[][4] = { {N3, N3, 0, 0}, /* 0 */ {P3, N3, 0, 0}, /* 1 */ {N3, P3, 0, 0}, /* 2 */ {P3, P3, 0, 0}, /* 3 */ {0, 0, 0, 0} /* backstop */ }; /* * Digits 0-2 (hour tens) */ #define P2 (P15 / 2) /* mark (+1) */ #define N2 (N15 / 2) /* space (-1) */ double bcd2[][4] = { {N2, N2, 0, 0}, /* 0 */ {P2, N2, 0, 0}, /* 1 */ {N2, P2, 0, 0}, /* 2 */ {0, 0, 0, 0} /* backstop */ }; /* * DST decode (DST2 DST1) for prettyprint */ char dstcod[] = { 'S', /* 00 standard time */ 'I', /* 01 daylight warning */ 'O', /* 10 standard warning */ '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 mldigit; /* maximum likelihood digit */ int phase; /* maximum likelihood digit phase */ int count; /* match count */ double digprb; /* max digit probability */ double digsnr; /* likelihood function (dB) */ double like[10]; /* likelihood integrator 0-9 */ }; /* * The station structure 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 amp; /* sync amplitude (I, Q square) */ double synamp; /* sync envelope at 800 ms */ double synmax; /* sync envelope at 0 s */ double synmin; /* avg sync envelope at 59 s, 1 s */ double synsnr; /* sync signal SNR */ double noise; /* max amplitude off pulse */ double sigmax; /* max amplitude on pulse */ double lastmax; /* last max amplitude on pulse */ long pos; /* position at maximum amplitude */ long lastpos; /* last position at maximum amplitude */ long jitter; /* shake, wiggle and waggle */ long mepoch; /* minute synch epoch */ int count; /* compare counter */ char refid[5]; /* reference identifier */ char ident[4]; /* station identifier */ int select; /* select bits */ }; /* * The channel structure is used to mitigate between channels. At this * point we have already decided which station to use. */ struct chan { int gain; /* audio gain */ int errcnt; /* data bit error counter */ double noiamp; /* I-channel average noise amplitude */ struct sync wwv; /* wwv station */ struct sync wwvh; /* wwvh station */ }; /* * WWV unit control structure */ struct wwvunit { l_fp timestamp; /* audio sample timestamp */ l_fp tick; /* audio sample increment */ double comp[SIZE]; /* decompanding table */ double phase, freq; /* logical clock phase and frequency */ double monitor; /* audio monitor point */ int fd_icom; /* ICOM file descriptor */ int errflg; /* error flags */ int bufcnt; /* samples in buffer */ int bufptr; /* buffer index pointer */ int port; /* codec port */ int gain; /* codec gain */ int clipcnt; /* sample clipped count */ int seccnt; /* second countdown */ int minset; /* minutes since last clock set */ int watch; /* watchcat */ int swatch; /* second sync watchcat */ /* * Variables used to establish basic system timing */ int avgint; /* log2 master time constant (s) */ int epoch; /* second epoch ramp */ int repoch; /* receiver sync epoch */ int yepoch; /* transmitter sync epoch */ double epomax; /* second sync amplitude */ 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; /* receiver sample counter */ int rsec; /* receiver seconds counter */ long mphase; /* minute sample counter */ long nepoch; /* minute epoch index */ /* * 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 cdelay; /* WWV propagation delay (samples) */ int hdelay; /* WVVH propagation delay (samples) */ int pdelay; /* propagation delay (samples) */ int tphase; /* transmitter sample counter */ int tsec; /* transmitter seconds counter */ int digcnt; /* count of digits synchronized */ /* * Variables used to estimate signal levels and bit/digit * probabilities */ double sigamp; /* I-channel peak signal amplitude */ double noiamp; /* I-channel average noise amplitude */ double datsnr; /* data 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 */ }; /* * Function prototypes */ static int wwv_start P((int, struct peer *)); static void wwv_shutdown P((int, struct peer *)); static void wwv_receive P((struct recvbuf *)); static void wwv_poll P((int, struct peer *)); /* * More function prototypes */ static void wwv_epoch P((struct peer *)); static void wwv_rf P((struct peer *, double)); static void wwv_endpoc P((struct peer *, double, int)); static void wwv_rsec P((struct peer *, double)); static void wwv_qrz P((struct peer *, struct sync *, double)); static void wwv_corr4 P((struct peer *, struct decvec *, double [], double [][4])); static void wwv_gain P((struct peer *)); static void wwv_tsec P((struct wwvunit *)); static double wwv_data P((struct wwvunit *, double)); static int timecode P((struct wwvunit *, char *)); static double wwv_snr P((double, double)); static int carry P((struct decvec *)); static void wwv_newchan P((struct peer *)); static int wwv_qsy P((struct peer *, int)); static double qsy[NCHAN] = {2.5, 5, 10, 15, 20}; /* frequencies (MHz) */ /* * 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 (not used) */ struct peer *peer /* peer structure pointer */ ) { struct refclockproc *pp; struct wwvunit *up; struct chan *cp; #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); if (fd < 0) return (0); #ifdef DEBUG if (debug) audio_show(); #endif /* * Allocate and initialize unit structure */ if (!(up = (struct wwvunit *) emalloc(sizeof(struct wwvunit)))) { (void) close(fd); return (0); } memset((char *)up, 0, sizeof(struct wwvunit)); pp = peer->procptr; pp->unitptr = (caddr_t)up; pp->io.clock_recv = wwv_receive; pp->io.srcclock = (caddr_t)peer; pp->io.datalen = 0; pp->io.fd = fd; if (!io_addclock(&pp->io)) { (void)close(fd); free(up); return (0); } /* * Initialize miscellaneous variables */ peer->precision = PRECISION; pp->clockdesc = DESCRIPTION; memcpy((char *)&pp->refid, REFID, 4); DTOLFP(1. / SECOND, &up->tick); /* * 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.; } /* * 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; /* * Initialize the station processes for audio gain, select bit, * station/frequency identifier and reference identifier. */ up->gain = 127; for (i = 0; i < NCHAN; i++) { cp = &up->mitig[i]; cp->gain = up->gain; cp->wwv.select = SELV; strcpy(cp->wwv.refid, "WWV "); sprintf(cp->wwv.ident,"C%.0f", floor(qsy[i])); cp->wwvh.select = SELH; strcpy(cp->wwvh.refid, "WWVH"); sprintf(cp->wwvh.ident, "H%.0f", floor(qsy[i])); } /* * Initialize autotune if available. Start out at 15 MHz. Note * that the ICOM select code must be less than 128, so the high * order bit can be used to select the line speed. */ #ifdef ICOM temp = 0; #ifdef DEBUG if (debug > 1) temp = P_TRACE; #endif if (peer->ttlmax != 0) { if (peer->ttlmax & 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) { up->schan = 3; if ((temp = wwv_qsy(peer, up->schan)) < 0) { NLOG(NLOG_SYNCEVENT | NLOG_SYSEVENT) msyslog(LOG_ERR, "ICOM bus error; autotune disabled"); up->errflg = CEVNT_FAULT; close(up->fd_icom); up->fd_icom = 0; } } #endif /* ICOM */ 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 = (struct wwvunit *)pp->unitptr; io_closeclock(&pp->io); if (up->fd_icom > 0) close(up->fd_icom); free(up); } /* * 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; /* * Local variables */ double sample; /* codec sample */ u_char *dpt; /* buffer pointer */ l_fp ltemp; int isneg; double dtemp; int i, j; peer = (struct peer *)rbufp->recv_srcclock; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; /* * Main loop - read until there ain't no more. Note codec * samples are bit-inverted. */ up->timestamp = rbufp->recv_time; up->bufcnt = rbufp->recv_length; DTOLFP((double)up->bufcnt / SECOND, <emp); L_SUB(&up->timestamp, <emp); dpt = rbufp->recv_buffer; for (up->bufptr = 0; up->bufptr < up->bufcnt; up->bufptr++) { sample = up->comp[~*dpt & 0xff]; /* * Clip noise spikes greater than MAXSIG. If no clips, * increase the gain a tad; if the clips are too high, * decrease a tad. */ if (sample > MAXSIG) { sample = MAXSIG; up->clipcnt++; } else if (sample < -MAXSIG) { sample = -MAXSIG; up->clipcnt++; } /* * Variable frequency oscillator. A phase change of one * unit produces a change of 360 degrees; a frequency * change of one unit produces a change of 1 Hz. */ up->phase += up->freq / SECOND; 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); /* * Once each second adjust the codec port and gain. * While at it, initialize the propagation delay for * both WWV and WWVH. Don't forget to correct for the * receiver phase delay, mostly due to the 600-Hz * IIR bandpass filter used for the sync signals. */ up->cdelay = (int)(SECOND * (pp->fudgetime1 + PDELAY)); up->hdelay = (int)(SECOND * (pp->fudgetime2 + PDELAY)); up->seccnt = (up->seccnt + 1) % SECOND; if (up->seccnt == 0) { if (pp->sloppyclockflag & CLK_FLAG2) up->port = 2; else up->port = 1; } /* * During development, it is handy to have an audio * monitor that can be switched to various signals. This * code converts the linear signal left in up->monitor * to codec format. */ isneg = 0; dtemp = up->monitor; if (sample < 0) { isneg = 1; dtemp -= dtemp; } i = 0; j = OFFSET >> 1; while (j != 0) { if (dtemp > up->comp[i]) i += j; else if (dtemp < up->comp[i]) i -= j; else break; j >>= 1; } if (isneg) *dpt = ~(i + OFFSET); else *dpt = ~i; dpt++; } /* * Squawk to the monitor speaker if enabled. */ if (pp->sloppyclockflag & CLK_FLAG3) if (write(pp->io.fd, (u_char *)&rbufp->recv_space, (u_int)up->bufcnt) < 0) perror("wwv:"); } /* * wwv_poll - called by the transmit procedure * * This routine keeps track of status. If nothing is heard for two * successive poll intervals, a timeout event is declared and any * orphaned timecode updates are sent to foster care. Once the clock is * set, it always appears reachable, unless reset by watchdog timeout. */ static void wwv_poll( int unit, /* instance number (not used) */ struct peer *peer /* peer structure pointer */ ) { struct refclockproc *pp; struct wwvunit *up; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; if (pp->coderecv == pp->codeproc) up->errflg = CEVNT_TIMEOUT; else pp->polls++; if (up->status & INSYNC) peer->reach |= 1; if (up->errflg) refclock_report(peer, up->errflg); up->errflg = 0; } /* * wwv_rf - process signals and demodulate to baseband * * This routine grooms and filters decompanded raw audio samples. The * output signals include the 100-Hz baseband data signal in quadrature * form, plus the epoch index of the second sync signal and the second * index of the minute sync signal. * * There are three 1-s ramps used by this program, all spanning the * range 0-7999 logical samples for exactly one second, as determined by * the logical clock. The first drives the second epoch and runs * continuously. The second determines the receiver phase and the third * the transmitter phase within the second. The receiver second begins * upon arrival of the 5-ms second sync pulse which begins the second; * while the transmitter second begins before it by the specified * propagation delay. * * There are three 1-m ramps spanning the range 0-59 seconds. The first * drives the minute epoch in samples and runs continuously. The second * determines the receiver second and the third the transmitter second. * The receiver second begins upon arrival of the 800-ms sync pulse sent * during the first second of the minute; while the transmitter second * begins before it by the specified propagation delay. * * The output signals include the epoch maximum and phase and second * maximum and index. The epoch phase provides the master reference for * all signal and timing functions, while the second index identifies * the first second of the minute. The epoch and second maxima are used * to calculate SNR for gating functions. * * Demodulation operations are based on three synthesized quadrature * sinusoids: 100 Hz for the data subcarrier, 1000 Hz for the WWV sync * signals and 1200 Hz for the WWVH sync signal. These drive synchronous * matched filters for the data subcarrier (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 seconds sync signal (5 ms at 1000 Hz) and * WWVH seconds sync signal (5 ms at 1200 Hz). */ static void wwv_rf( struct peer *peer, /* peerstructure pointer */ double isig /* input signal */ ) { struct refclockproc *pp; struct wwvunit *up; static double lpf[5]; /* 150-Hz lpf delay line */ double data; /* lpf output */ static double bpf[9]; /* 1000/1200-Hz bpf delay line */ double syncx; /* bpf output */ static double mf[41]; /* 1000/1200-Hz mf delay line */ double mfsync; /* mf output */ static int iptr; /* data channel pointer */ static double ibuf[DATSIZ]; /* data I channel delay line */ static double qbuf[DATSIZ]; /* data Q channel delay line */ static int jptr; /* sync channel pointer */ static 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 int csinptr; /* wwv 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 int hsinptr; /* wwvh channels phase */ static double epobuf[SECOND]; /* epoch sync comb filter */ static double epomax; /* epoch sync amplitude buffer */ static int epopos; /* epoch sync position buffer */ static int iniflg; /* initialization flag */ struct sync *sp; double dtemp; long ltemp; int i; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; if (!iniflg) { iniflg = 1; memset((char *)lpf, 0, sizeof(lpf)); memset((char *)bpf, 0, sizeof(bpf)); memset((char *)mf, 0, sizeof(mf)); memset((char *)ibuf, 0, sizeof(ibuf)); memset((char *)qbuf, 0, sizeof(qbuf)); memset((char *)cibuf, 0, sizeof(cibuf)); memset((char *)cqbuf, 0, sizeof(cqbuf)); memset((char *)hibuf, 0, sizeof(hibuf)); memset((char *)hqbuf, 0, sizeof(hqbuf)); memset((char *)epobuf, 0, sizeof(epobuf)); } up->monitor = isig; /* change for debug */ /* * 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. * * Matlab IIR 4th-order IIR elliptic, 150 Hz lowpass, 0.2 dB * passband ripple, -50 dB stopband ripple. */ data = (lpf[4] = lpf[3]) * 8.360961e-01; data += (lpf[3] = lpf[2]) * -3.481740e+00; data += (lpf[2] = lpf[1]) * 5.452988e+00; data += (lpf[1] = lpf[0]) * -3.807229e+00; lpf[0] = isig - 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 I and Q quadrature data signals are produced by * multiplying the filtered signal by 100-Hz sine and cosine * signals, respectively. The data signals are demodulated by * 170-ms synchronous matched filters to produce the amplitude * and phase signals used by the decoder. Note the correction * due to the propagation delay is necessary for seamless * handover between WWV and WWVH. */ i = up->datapt - up->pdelay % 80; if (i < 0) i += 80; up->datapt = (up->datapt + IN100) % 80; dtemp = sintab[i] * data / DATSIZ * DGAIN; up->irig -= ibuf[iptr]; ibuf[iptr] = dtemp; up->irig += dtemp; i = (i + 20) % 80; dtemp = sintab[i] * data / DATSIZ * DGAIN; up->qrig -= qbuf[iptr]; qbuf[iptr] = dtemp; up->qrig += dtemp; iptr = (iptr + 1) % DATSIZ; /* * Baseband sync demodulation. The 1000/1200 sync signals are * extracted using a 600-Hz IIR bandpass filter. This removes * the 100-Hz data subcarrier, as well as the 440-Hz and 600-Hz * tones and most of the noise and voice modulation components. * * Matlab 4th-order IIR elliptic, 800-1400 Hz bandpass, 0.2 dB * passband ripple, -50 dB stopband ripple. */ syncx = (bpf[8] = bpf[7]) * 4.897278e-01; syncx += (bpf[7] = bpf[6]) * -2.765914e+00; syncx += (bpf[6] = bpf[5]) * 8.110921e+00; syncx += (bpf[5] = bpf[4]) * -1.517732e+01; syncx += (bpf[4] = bpf[3]) * 1.975197e+01; syncx += (bpf[3] = bpf[2]) * -1.814365e+01; syncx += (bpf[2] = bpf[1]) * 1.159783e+01; syncx += (bpf[1] = bpf[0]) * -4.735040e+00; bpf[0] = isig - syncx; syncx = bpf[0] * 8.203628e-03 + bpf[1] * -2.375732e-02 + bpf[2] * 3.353214e-02 + bpf[3] * -4.080258e-02 + bpf[4] * 4.605479e-02 + bpf[5] * -4.080258e-02 + bpf[6] * 3.353214e-02 + bpf[7] * -2.375732e-02 + bpf[8] * 8.203628e-03; /* * The I and Q quadrature minute sync signals are produced by * multiplying the filtered signal by 1000-Hz (WWV) and 1200-Hz * (WWVH) sine and cosine signals, respectively. The resulting * signals are demodulated by 800-ms synchronous matched filters * to synchronize the second and minute and to detect which one * (or both) the WWV or WWVH signal is present. */ up->mphase = (up->mphase + 1) % MINUTE; i = csinptr; csinptr = (csinptr + IN1000) % 80; dtemp = sintab[i] * syncx / SYNSIZ * SGAIN; ciamp = ciamp - cibuf[jptr] + dtemp; cibuf[jptr] = dtemp; i = (i + 20) % 80; dtemp = sintab[i] * syncx / SYNSIZ * SGAIN; cqamp = cqamp - cqbuf[jptr] + dtemp; cqbuf[jptr] = dtemp; dtemp = ciamp * ciamp + cqamp * cqamp; wwv_qrz(peer, &up->mitig[up->schan].wwv, dtemp); i = hsinptr; hsinptr = (hsinptr + IN1200) % 80; dtemp = sintab[i] * syncx / SYNSIZ * SGAIN; hiamp = hiamp - hibuf[jptr] + dtemp; hibuf[jptr] = dtemp; i = (i + 20) % 80; dtemp = sintab[i] * syncx / SYNSIZ * SGAIN; hqamp = hqamp - hqbuf[jptr] + dtemp; hqbuf[jptr] = dtemp; dtemp = hiamp * hiamp + hqamp * hqamp; wwv_qrz(peer, &up->mitig[up->schan].wwvh, dtemp); jptr = (jptr + 1) % SYNSIZ; if (up->mphase == 0) { /* * This section is called once per minute at the minute * epoch independently of the transmitter or receiver * minute. If the leap bit is set, set the minute epoch * back one second so the station processes don't miss a * beat. Then, increment the watchdog counter and test * for two sets of conditions depending on whether * minute sync has been acquired or not. */ up->watch++; if (up->rsec == 60) { up->mphase -= SECOND; if (up->mphase < 0) up->mphase += MINUTE; } else if (!(up->status & MSYNC)) { /* * If minute sync has not been acquired, the * program listens for minute sync pulses from * both WWV and WWVH. The station with the * greater compare count is selected, with ties * broken by WWV, but only if the count is at * least three. Once a station has been * acquired, it is initialized and begins * tracking the signal. */ if (up->mitig[up->achan].wwv.count >= up->mitig[up->achan].wwvh.count) sp = &up->mitig[up->achan].wwv; else sp = &up->mitig[up->achan].wwvh; if (sp->count >= AMIN) { up->watch = up->swatch = 0; up->status |= MSYNC; ltemp = sp->mepoch - SYNSIZ; if (ltemp < 0) ltemp += MINUTE; up->rsec = (MINUTE - ltemp) / SECOND; if (!(up->status & SSYNC)) { up->repoch = ltemp % SECOND; up->yepoch = up->repoch - up->pdelay; if (up->yepoch < 0) up->yepoch += SECOND; } wwv_newchan(peer); } else if (sp->count == 0 || up->watch >= ACQSN) { up->watch = sp->count = 0; up->schan = (up->schan + 1) % NCHAN; wwv_qsy(peer, up->schan); } } else { /* * If minute sync has been acquired, the program * watches for timeout. The timeout is reset * when the clock is set or verified. If a * timeout occurs and the minute units digit has * not synchronized, reset the program and start * over. */ if (up->watch > DIGIT && !(up->status & DSYNC)) up->watch = up->status = 0; /* * If the second sync times out, dim the sync * lamp and raise an alarm. */ up->swatch++; if (up->swatch > HSPEC) up->status &= ~SSYNC; if (!(up->status & SSYNC)) up->alarm |= 1 << SYNERR; } } /* * The second sync pulse is extracted using 5-ms 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). */ if (up->status & SELV) { up->pdelay = up->cdelay; /* * WWV FIR matched filter, five cycles of 1000-Hz * sinewave. */ mf[40] = mf[39]; mfsync = (mf[39] = mf[38]) * 4.224514e-02; mfsync += (mf[38] = mf[37]) * 5.974365e-02; mfsync += (mf[37] = mf[36]) * 4.224514e-02; mf[36] = mf[35]; mfsync += (mf[35] = mf[34]) * -4.224514e-02; mfsync += (mf[34] = mf[33]) * -5.974365e-02; mfsync += (mf[33] = mf[32]) * -4.224514e-02; mf[32] = mf[31]; mfsync += (mf[31] = mf[30]) * 4.224514e-02; mfsync += (mf[30] = mf[29]) * 5.974365e-02; mfsync += (mf[29] = mf[28]) * 4.224514e-02; mf[28] = mf[27]; mfsync += (mf[27] = mf[26]) * -4.224514e-02; mfsync += (mf[26] = mf[25]) * -5.974365e-02; mfsync += (mf[25] = mf[24]) * -4.224514e-02; mf[24] = mf[23]; mfsync += (mf[23] = mf[22]) * 4.224514e-02; mfsync += (mf[22] = mf[21]) * 5.974365e-02; mfsync += (mf[21] = mf[20]) * 4.224514e-02; mf[20] = mf[19]; mfsync += (mf[19] = mf[18]) * -4.224514e-02; mfsync += (mf[18] = mf[17]) * -5.974365e-02; mfsync += (mf[17] = mf[16]) * -4.224514e-02; mf[16] = mf[15]; mfsync += (mf[15] = mf[14]) * 4.224514e-02; mfsync += (mf[14] = mf[13]) * 5.974365e-02; mfsync += (mf[13] = mf[12]) * 4.224514e-02; mf[12] = mf[11]; mfsync += (mf[11] = mf[10]) * -4.224514e-02; mfsync += (mf[10] = mf[9]) * -5.974365e-02; mfsync += (mf[9] = mf[8]) * -4.224514e-02; mf[8] = mf[7]; mfsync += (mf[7] = mf[6]) * 4.224514e-02; mfsync += (mf[6] = mf[5]) * 5.974365e-02; mfsync += (mf[5] = mf[4]) * 4.224514e-02; mf[4] = mf[3]; mfsync += (mf[3] = mf[2]) * -4.224514e-02; mfsync += (mf[2] = mf[1]) * -5.974365e-02; mfsync += (mf[1] = mf[0]) * -4.224514e-02; mf[0] = syncx; } else if (up->status & SELH) { up->pdelay = up->hdelay; /* * WWVH FIR matched filter, six cycles of 1200-Hz * sinewave. */ mf[40] = mf[39]; mfsync = (mf[39] = mf[38]) * 4.833363e-02; mfsync += (mf[38] = mf[37]) * 5.681959e-02; mfsync += (mf[37] = mf[36]) * 1.846180e-02; mfsync += (mf[36] = mf[35]) * -3.511644e-02; mfsync += (mf[35] = mf[34]) * -5.974365e-02; mfsync += (mf[34] = mf[33]) * -3.511644e-02; mfsync += (mf[33] = mf[32]) * 1.846180e-02; mfsync += (mf[32] = mf[31]) * 5.681959e-02; mfsync += (mf[31] = mf[30]) * 4.833363e-02; mf[30] = mf[29]; mfsync += (mf[29] = mf[28]) * -4.833363e-02; mfsync += (mf[28] = mf[27]) * -5.681959e-02; mfsync += (mf[27] = mf[26]) * -1.846180e-02; mfsync += (mf[26] = mf[25]) * 3.511644e-02; mfsync += (mf[25] = mf[24]) * 5.974365e-02; mfsync += (mf[24] = mf[23]) * 3.511644e-02; mfsync += (mf[23] = mf[22]) * -1.846180e-02; mfsync += (mf[22] = mf[21]) * -5.681959e-02; mfsync += (mf[21] = mf[20]) * -4.833363e-02; mf[20] = mf[19]; mfsync += (mf[19] = mf[18]) * 4.833363e-02; mfsync += (mf[18] = mf[17]) * 5.681959e-02; mfsync += (mf[17] = mf[16]) * 1.846180e-02; mfsync += (mf[16] = mf[15]) * -3.511644e-02; mfsync += (mf[15] = mf[14]) * -5.974365e-02; mfsync += (mf[14] = mf[13]) * -3.511644e-02; mfsync += (mf[13] = mf[12]) * 1.846180e-02; mfsync += (mf[12] = mf[11]) * 5.681959e-02; mfsync += (mf[11] = mf[10]) * 4.833363e-02; mf[10] = mf[9]; mfsync += (mf[9] = mf[8]) * -4.833363e-02; mfsync += (mf[8] = mf[7]) * -5.681959e-02; mfsync += (mf[7] = mf[6]) * -1.846180e-02; mfsync += (mf[6] = mf[5]) * 3.511644e-02; mfsync += (mf[5] = mf[4]) * 5.974365e-02; mfsync += (mf[4] = mf[3]) * 3.511644e-02; mfsync += (mf[3] = mf[2]) * -1.846180e-02; mfsync += (mf[2] = mf[1]) * -5.681959e-02; mfsync += (mf[1] = mf[0]) * -4.833363e-02; mf[0] = syncx; } else { mfsync = 0; } /* * Extract the seconds sync pulse using a 1-s comb filter at * baseband. Correct for the FIR matched filter delay, which is * 5 ms for both the WWV and WWVH filters. Blank the signal when * probing. */ up->epoch = (up->epoch + 1) % SECOND; if (up->epoch == 0) { wwv_endpoc(peer, epomax, epopos); up->epomax = epomax; epomax = 0; if (!(up->status & MSYNC)) wwv_gain(peer); } dtemp = (epobuf[up->epoch] += (mfsync - epobuf[up->epoch]) / (MINAVG << up->avgint)); if (dtemp > epomax) { epomax = dtemp; epopos = up->epoch - up->pdelay - 5 * MS; if (epopos < 0) epopos += SECOND; } if (up->status & MSYNC) wwv_epoch(peer); } /* * 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 in turn for the minute pulse from either station, which * involves searching through the entire epoch minute of samples. After * minute sync acquisition, the process searches only during the probe * window, which occupies seconds 59, 0 and 1, to construct a metric * used to determine which frequency and station provides the best * signal. * * The pulse discriminator requires that (a) the peak on-pulse sample * amplitude must be above 2000, (b) the SNR relative to the peak * off-pulse sample amplitude must be reduced 6 dB or more below the * peak and (c) the maximum difference between the current and previous * epoch indices must be less than 50 ms. A compare counter keeps track * of the number of successive intervals which satisfy these criteria. * * Students of radar receiver technology will discover this algorithm * amounts to a range gate discriminator. In practice, the performance * of this gadget is amazing. Once setting teeth in a station, it hangs * on until the minute beep can barely be heard and long after the * second tick and comb filter have given up. */ static void wwv_qrz( struct peer *peer, /* peerstructure pointer */ struct sync *sp, /* sync channel structure */ double syncx /* bandpass filtered sync signal */ ) { struct refclockproc *pp; struct wwvunit *up; char tbuf[80]; /* monitor buffer */ double snr; /* on-pulse/off-pulse ratio (dB) */ long epoch; int isgood; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; /* * Find the sample with peak energy, which defines the minute * epoch. If minute sync has been acquired, search only the * probe window; otherwise, search the entire minute. If a * maximum has been found with good amplitude, search only the * second before and after that position for the next maximum * and the rest of the window for the noise. */ if (!(up->status & MSYNC) || up->status & SFLAG) { sp->amp = syncx; if (up->status & MSYNC) epoch = up->nepoch; else if (sp->count > 1) epoch = sp->mepoch; else epoch = sp->lastpos; if (syncx > sp->sigmax) { sp->sigmax = syncx; sp->pos = up->mphase; } if (abs(MOD(up->mphase - epoch, MINUTE)) > SYNSIZ && syncx > sp->noise) { sp->noise = syncx; } } if (up->mphase == 0) { /* * At the end of the minute, determine the epoch of the * sync pulse, as well as the SNR and difference between * the current and previous epoch (jitter). */ sp->jitter = MOD(sp->pos - sp->lastpos, MINUTE); sp->select &= ~JITRNG; if (abs(sp->jitter) > AWND * MS) sp->select |= JITRNG; sp->sigmax = SQRT(sp->sigmax); sp->noise = SQRT(sp->noise); if (up->status & MSYNC) { /* * If in minute sync, just count the runs up and * down. */ if (sp->select & (DATANG | SYNCNG | JITRNG)) { if (sp->count > 0) sp->count--; } else { if (sp->count < AMAX) sp->count++; } } else { /* * If not yet in minute sync, we have to do a * little dance to find a valid minute sync * pulse, emphasis valid. */ snr = wwv_snr(sp->sigmax, sp->noise); isgood = sp->sigmax > ATHR && snr > ASNR && !(sp->select & JITRNG); switch (sp->count) { /* * In state 0 the station was not heard during * the previous probe. Look for the biggest blip * greater than the amplitude threshold in the * minute and assume that the minute sync pulse. * If found, bump to state 1. */ case 0: if (sp->sigmax >= ATHR) sp->count++; break; /* * In state 1 a candidate blip has been found * and the next minute has been searched for * another blip. If none are found greater than * the threshold, or if the biggest blip outside * the candidate pulse is less than 6 dB below * the biggest blip, drop back to state 0 and * hunt some more. Otherwise, a legitimate * minute pulse may have been found, so bump to * state 2. */ case 1: if (sp->sigmax < ATHR) { sp->count--; break; } else if (!isgood) { break; } /* fall through */ /* * In states 2 and above, continue to groom * samples as before and drop back to the * previous state if the groom fails. If it * succeeds, bump to the next state until * reaching the clamp, if ever. */ default: if (!isgood) { sp->count--; break; } sp->mepoch = sp->pos; if (sp->count < AMAX) sp->count++; break; } sprintf(tbuf, "wwv8 %d %3d %-3s %d %5.0f %5.1f %7ld %7ld %7ld", up->port, up->gain, sp->ident, sp->count, sp->sigmax, snr, sp->pos, sp->jitter, MOD(sp->pos - up->nepoch - SYNSIZ, MINUTE)); if (pp->sloppyclockflag & CLK_FLAG4) record_clock_stats(&peer->srcadr, tbuf); #ifdef DEBUG if (debug) printf("%s\n", tbuf); #endif } sp->lastmax = sp->sigmax; sp->lastpos = sp->pos; sp->sigmax = sp->noise = 0; } } /* * wwv_endpoc - process receiver epoch * * This routine is called at the end of the receiver epoch. It * determines the epoch position within the second and disciplines the * sample clock using a frequency-lock loop (FLL). * * Seconds 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 data filtering and grooming. */ static void wwv_endpoc( struct peer *peer, /* peer structure pointer */ double epomax, /* epoch max */ int epopos /* epoch max position */ ) { struct refclockproc *pp; struct wwvunit *up; static int epoch_mf[3]; /* epoch median filter */ static int tepoch; /* median filter epoch */ static int tspan; /* median filter span */ static int xepoch; /* last second epoch */ static int zepoch; /* last averaging interval epoch */ static int syncnt; /* second epoch run length counter */ static int jitcnt; /* jitter holdoff counter */ static int avgcnt; /* averaging interval counter */ static int avginc; /* averaging ratchet */ static int iniflg; /* initialization flag */ char tbuf[80]; /* monitor buffer */ double dtemp; int tmp2, tmp3; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; if (!iniflg) { iniflg = 1; memset((char *)epoch_mf, 0, sizeof(epoch_mf)); } /* * A three-stage median filter is used to help denoise the * seconds sync pulse. The median sample becomes the candidate * epoch; the difference between the other two samples becomes * the span, which is used currently only for debugging. */ 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 */ tspan = epoch_mf[0] - epoch_mf[2]; } else if (epoch_mf[2] > epoch_mf[0]) { tepoch = epoch_mf[0]; /* 2 0 1 */ tspan = epoch_mf[2] - epoch_mf[1]; } else { tepoch = epoch_mf[2]; /* 0 2 1 */ tspan = epoch_mf[0] - epoch_mf[1]; } } else { if (epoch_mf[1] < epoch_mf[2]) { tepoch = epoch_mf[1]; /* 2 1 0 */ tspan = epoch_mf[2] - epoch_mf[0]; } else if (epoch_mf[2] < epoch_mf[0]) { tepoch = epoch_mf[0]; /* 1 0 2 */ tspan = epoch_mf[1] - epoch_mf[2]; } else { tepoch = epoch_mf[2]; /* 1 2 0 */ tspan = epoch_mf[1] - epoch_mf[0]; } } /* * If the epoch candidate is within 1 ms of the last one, the * new candidate replaces the last one and the jitter counter is * reset; otherwise, the candidate is ignored and the jitter * counter is incremented. If the jitter counter exceeds the * frequency averaging interval, the new candidate replaces the * old one anyway. The compare counter is incremented if the new * candidate is identical to the last one; otherwise, it is * forced to zero. If the compare counter increments to 10, the * epoch is reset and the receiver second epoch is set. * * Careful attention to detail here. If the signal amplitude * falls below the threshold or if no stations are heard, we * certainly cannot be in sync. */ tmp2 = MOD(tepoch - xepoch, SECOND); if (up->epomax < STHR || !(up->status & (SELV | SELH))) { up->status &= ~SSYNC; jitcnt = syncnt = avgcnt = 0; } else if (abs(tmp2) <= MS || jitcnt >= (MINAVG << up->avgint)) { jitcnt = 0; if (tmp2 != 0) { xepoch = tepoch; syncnt = 0; } else { if (syncnt < SCMP) { syncnt++; } else { up->status |= SSYNC; up->swatch = 0; up->repoch = tepoch; up->yepoch = up->repoch; if (up->yepoch < 0) up->yepoch += SECOND; } } avgcnt++; } else { jitcnt++; syncnt = avgcnt = 0; } if (!(up->status & SSYNC) && 0) { sprintf(tbuf, "wwv1 %2d %04x %5.0f %2d %5.0f %5d %5d %5d %2d %4d", up->rsec, up->status, up->epomax, avgcnt, epomax, tepoch, tspan, tmp2, syncnt, jitcnt); if (pp->sloppyclockflag & CLK_FLAG4) record_clock_stats(&peer->srcadr, tbuf); #ifdef DEBUG if (debug) printf("%s\n", tbuf); #endif /* DEBUG */ } /* * The sample clock frequency is disciplined using a first-order * feedback loop with time constant consistent with the Allan * intercept of typical computer clocks. The loop update is * calculated each averaging interval from the epoch change in * 125-us units and interval length in seconds. The interval is * doubled after four intervals where epoch change is not more * than one sample. * * The averaging interval affects other receiver functions, * including the the 1000/1200-Hz comb filter and sample clock * loop. It also affects the 100-Hz subcarrier loop and the bit * and digit comparison counter thresholds. */ tmp3 = MOD(tepoch - zepoch, SECOND); if (avgcnt >= (MINAVG << up->avgint)) { if (abs(tmp3) < MS) { dtemp = (double)tmp3 / avgcnt; up->freq += dtemp / SYNCTC; if (up->freq > MAXFREQ) up->freq = MAXFREQ; else if (up->freq < -MAXFREQ) up->freq = -MAXFREQ; if (abs(tmp3) <= 1 && up->avgint < MAXAVG) { if (avginc < 4) { avginc++; } else { avginc = 0; up->avgint++; } } if (up->avgint < MAXAVG) { sprintf(tbuf, "wwv2 %2d %04x %5.0f %5d %5d %2d %2d %6.1f %6.1f", up->rsec, up->status, up->epomax, MINAVG << up->avgint, avgcnt, avginc, tmp3, dtemp / SECOND * 1e6, up->freq / SECOND * 1e6); if (pp->sloppyclockflag & CLK_FLAG4) record_clock_stats( &peer->srcadr, tbuf); #ifdef DEBUG if (debug) printf("%s\n", tbuf); #endif /* DEBUG */ } } zepoch = tepoch; avgcnt = 0; } } /* * wwv_epoch - main loop * * This routine establishes receiver and transmitter epoch * synchronization and determines the data subcarrier pulse length. * 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. This establishes when to sample the data * subcarrier in-phase signal for the maximum level and noise level and * when to determine the pulse length. The transmitter second leads the * receiver second by the propagation delay, receiver delay and filter * delay of this program. It establishes the clock time and implements * the sometimes idiosyncratic conventional clock time and civil * calendar. * * 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 */ ) { static double dpulse; /* data pulse length */ struct refclockproc *pp; struct wwvunit *up; struct chan *cp; struct sync *sp; l_fp offset; /* NTP format offset */ double dtemp; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; /* * Sample the minute sync pulse amplitude at epoch 800 for both * the WWV and WWVH stations. This will be used later for * channel mitigation. */ cp = &up->mitig[up->achan]; if (up->rphase == 800 * MS) { sp = &cp->wwv; sp->synamp = SQRT(sp->amp); sp = &cp->wwvh; sp->synamp = SQRT(sp->amp); } if (up->rsec == 0) { up->sigamp = up->datsnr = 0; } else { /* * Estimate the noise level by integrating the I-channel * energy at epoch 30 ms. */ if (up->rphase == 30 * MS) { if (!(up->status & SFLAG)) up->noiamp += (up->irig - up->noiamp) / (MINAVG << up->avgint); else cp->noiamp += (SQRT(up->irig * up->irig + up->qrig * up->qrig) - cp->noiamp) / 8; /* * Strobe the peak I-channel data signal at epoch 200 * ms. Compute the SNR and adjust the 100-Hz reference * oscillator phase using the Q-channel data signal at * that epoch. Save the envelope amplitude for the probe * channel. */ } else if (up->rphase == 200 * MS) { if (!(up->status & SFLAG)) { up->sigamp = up->irig; if (up->sigamp < 0) up->sigamp = 0; up->datsnr = wwv_snr(up->sigamp, up->noiamp); up->datpha = up->qrig / (MINAVG << 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; } } else { up->sigamp = SQRT(up->irig * up->irig + up->qrig * up->qrig); up->datsnr = wwv_snr(up->sigamp, cp->noiamp); } /* * The slice level is set half way between the peak * signal and noise levels. Strobe the negative zero * crossing after epoch 200 ms and record the epoch at * that time. This defines the length of the data pulse, * which will later be converted into scaled bit * probabilities. */ } else if (up->rphase > 200 * MS) { dtemp = (up->sigamp + up->noiamp) / 2; if (up->irig < dtemp && dpulse == 0) dpulse = up->rphase; } } /* * At the end of the transmitter second, crank the clock state * machine. Note we have to be careful to set the transmitter * epoch at the same time as the receiver epoch to be sure the * right propagation delay is used. We don't bother the heavy * machinery unless the clock is set. */ up->tphase++; if (up->epoch == up->yepoch) { wwv_tsec(up); up->tphase = 0; /* * Determine the current offset from the time of century * and the sample timestamp, but only if the SYNERR * alarm has not been raised in the present or previous * minute. */ if (!(up->status & SFLAG) && up->status & INSYNC && (up->alarm & (3 << SYNERR)) == 0) { pp->second = up->tsec; pp->minute = up->decvec[MN].digit + up->decvec[MN + 1].digit * 10; pp->hour = up->decvec[HR].digit + up->decvec[HR + 1].digit * 10; pp->day = up->decvec[DA].digit + up->decvec[DA + 1].digit * 10 + up->decvec[DA + 2].digit * 100; pp->year = up->decvec[YR].digit + up->decvec[YR + 1].digit * 10; if (pp->year < UTCYEAR) pp->year += 2000; else pp->year += 1900; /* * We have to simulate refclock_process() here, * since the fudgetime gets added much earlier * than this. */ pp->lastrec = up->timestamp; L_CLR(&offset); if (!clocktime(pp->day, pp->hour, pp->minute, pp->second, GMT, pp->lastrec.l_ui, &pp->yearstart, &offset.l_ui)) up->errflg = CEVNT_BADTIME; else refclock_process_offset(pp, offset, pp->lastrec, 0.); } } /* * At the end of the receiver second, process the data bit and * update the decoding matrix probabilities. */ up->rphase++; if (up->epoch == up->repoch) { wwv_rsec(peer, dpulse); wwv_gain(peer); up->rphase = dpulse = 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. BTW, stations WWV/WWVH cowardly kill the * transmitter carrier for a few seconds around the leap to avoid icky * details of transmission format during the leap. */ static void wwv_rsec( struct peer *peer, /* peer structure pointer */ double dpulse ) { 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; double bit; /* bit likelihood */ char tbuf[80]; /* monitor buffer */ int sw, arg, nsec; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; if (!iniflg) { iniflg = 1; memset((char *)bitvec, 0, sizeof(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 + 1; bit = wwv_data(up, dpulse); bitvec[up->rsec] += (bit - bitvec[up->rsec]) / TCONST; sw = progx[up->rsec].sw; arg = progx[up->rsec].arg; 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 (not used), and the minute sync * pulse in second 0. At the end of second 58, we QSYed to the * probe channel, which rotates over all WWV/H frequencies. At * the end of second 59, we latched the sync noise and tested * for data bit error. At the end of second 0, we now latch the * sync peak. */ case SYNC2: /* 0 */ cp = &up->mitig[up->achan]; sp = &cp->wwv; sp->synmax = sp->synamp; sp = &cp->wwvh; sp->synmax = sp->synamp; break; /* * At the end of second 1, latch and average the sync noise and * test for data bit error. Set SYNCNG if the sync pulse * amplitude and SNR are not above thresholds. Set DATANG if * data error occured on both second 59 and second 1. Finally, * QSY back to the data channel. */ case SYNC3: /* 1 */ cp = &up->mitig[up->achan]; if (up->sigamp < DTHR || up->datsnr < DSNR) cp->errcnt++; sp = &cp->wwv; sp->synmin = (sp->synmin + sp->synamp) / 2; sp->synsnr = wwv_snr(sp->synmax, sp->synmin); sp->select &= ~(DATANG | SYNCNG); if (sp->synmax < QTHR || sp->synsnr < QSNR) sp->select |= SYNCNG; if (cp->errcnt > 1) sp->select |= DATANG; rp = &cp->wwvh; rp->synmin = (rp->synmin + rp->synamp) / 2; rp->synsnr = wwv_snr(rp->synmax, rp->synmin); rp->select &= ~(DATANG | SYNCNG); if (rp->synmax < QTHR || rp->synsnr < QSNR) rp->select |= SYNCNG; if (cp->errcnt > 1) rp->select |= DATANG; cp->errcnt = 0; sprintf(tbuf, "wwv5 %d %3d %-3s %04x %d %.0f/%.1f/%ld %s %04x %d %.0f/%.1f/%ld", up->port, up->gain, sp->ident, sp->select, sp->count, sp->synmax, sp->synsnr, sp->jitter, rp->ident, rp->select, rp->count, rp->synmax, rp->synsnr, rp->jitter); if (pp->sloppyclockflag & CLK_FLAG4) record_clock_stats(&peer->srcadr, tbuf); #ifdef DEBUG if (debug) printf("%s\n", tbuf); #endif /* DEBUG */ up->status &= ~SFLAG; wwv_newchan(peer); break; /* * Save the bit probability in the BCD data vector at the index * given by the argument. Note that all bits of the vector have * to be above the data gate threshold for the digit to be * considered valid. Bits not used in the digit are forced to * zero and not checked for errors. */ case COEF1: /* 10-13 */ if (up->status & DGATE) up->status |= BGATE; bcddld[arg] = bit; break; case COEF2: /* 18, 27-28, 42-43 */ bcddld[arg] = 0; break; case COEF: /* 4-7, 15-17, 20-23, 25-26, 30-33, 35-38, 40-41, 51-54 */ if (up->status & DGATE || !(up->status & DSYNC)) up->status |= BGATE; 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; case DECIM9: /* 8, 14, 24, 34, 39 */ wwv_corr4(peer, &up->decvec[arg], bcddld, bcd9); break; /* * Miscellaneous bits. If above the positive threshold, declare * 1; if below the negative threshold, declare 0; otherwise * raise the SYMERR alarm. At the end of second 58, QSY to the * probe channel. */ case MSC20: /* 55 */ wwv_corr4(peer, &up->decvec[YR + 1], bcddld, bcd9); /* fall through */ case MSCBIT: /* 2, 3, 50, 56-57 */ if (bitvec[up->rsec] > BTHR) up->misc |= arg; else if (bitvec[up->rsec] < -BTHR) up->misc &= ~arg; else up->alarm |= 1 << SYMERR; break; case MSC21: /* 58 */ if (bitvec[up->rsec] > BTHR) up->misc |= arg; else if (bitvec[up->rsec] < -BTHR) up->misc &= ~arg; else up->alarm |= 1 << SYMERR; up->schan = (up->schan + 1) % NCHAN; wwv_qsy(peer, up->schan); up->status |= SFLAG; break; /* * The endgames * * Second 59 contains the first data pulse of the probe * sequence. Check it for validity and establish the noise floor * for the minute sync SNR. */ case MIN1: /* 59 */ cp = &up->mitig[up->achan]; if (up->sigamp < DTHR || up->datsnr < DSNR) cp->errcnt++; sp = &cp->wwv; sp->synmin = sp->synamp; sp = &cp->wwvh; sp->synmin = sp->synamp; /* * If SECWARN is set on the last minute of 30 June or 31 * December, LEPSEC bit is set. At the end of the minute * in which LEPSEC is set the transmitter and receiver * insert an extra second (60) in the timescale and the * minute sync skips a second. We only get to test this * wrinkle at intervals of about 18 months, the actual * mileage may vary. */ if (up->tsec == 60) { up->status &= ~LEPSEC; break; } /* fall through */ /* * If all nine clock digits are valid and the SYNERR alarm is * not raised in the current or previous second, the clock is * set or validated. If at least one digit is set, which by * design must be the minute units digit, the clock state * machine begins to count the minutes. */ case MIN2: /* 59/60 */ up->minset = ((current_time - peer->update) + 30) / 60; if (up->digcnt > 0) up->status |= DSYNC; if (up->digcnt >= 9 && (up->alarm & (3 << SYNERR)) == 0) { up->status |= INSYNC; up->watch = 0; } pp->lencode = timecode(up, pp->a_lastcode); if (up->misc & SECWAR) pp->leap = LEAP_ADDSECOND; else pp->leap = LEAP_NOWARNING; refclock_receive(peer); record_clock_stats(&peer->srcadr, pp->a_lastcode); #ifdef DEBUG if (debug) printf("wwv: timecode %d %s\n", pp->lencode, pp->a_lastcode); #endif /* DEBUG */ /* * The ultimate watchdog is the interval since the * reference clock interface code last received an * update from this driver. If the interval is greater * than a couple of days, manual intervention is * probably required, so the program resets and tries to * resynchronized from scratch. */ if (up->minset > PANIC) up->status = 0; up->alarm = (up->alarm & ~0x8888) << 1; up->nepoch = (up->mphase + SYNSIZ) % MINUTE; up->errcnt = up->digcnt = nsec = 0; break; } if (!(up->status & DSYNC)) { sprintf(tbuf, "wwv3 %2d %04x %5.0f %5.0f %5.0f %5.1f %5.0f %5.0f", up->rsec, up->status, up->epomax, up->sigamp, up->datpha, up->datsnr, bit, bitvec[up->rsec]); if (pp->sloppyclockflag & CLK_FLAG4) record_clock_stats(&peer->srcadr, tbuf); #ifdef DEBUG if (debug) printf("%s\n", tbuf); #endif /* DEBUG */ } up->rsec = up->tsec = nsec; return; } /* * wwv_data - calculate bit probability * * This routine is called at the end of the receiver second to calculate * the probabilities that the previous second contained a zero (P0), one * (P1) or position indicator (P2) pulse. If not in sync or if the data * bit is bad, a bit error is declared and the probabilities are forced * to zero. Otherwise, the data gate is enabled and the probabilities * are calculated. Note that the data matched filter contributes half * the pulse width, or 85 ms.. */ static double wwv_data( struct wwvunit *up, /* driver unit pointer */ double pulse /* pulse length (sample units) */ ) { double p0, p1, p2; /* probabilities */ double dpulse; /* pulse length in ms */ p0 = p1 = p2 = 0; dpulse = pulse - DATSIZ / 2; /* * If the data amplitude or SNR are below threshold or if the * pulse length is less than 170 ms, declare an erasure. */ if (up->sigamp < DTHR || up->datsnr < DSNR || dpulse < DATSIZ) { up->status |= DGATE; up->errcnt++; if (up->errcnt > MAXERR) up->alarm |= 1 << MODERR; return (0); } /* * The probability of P0 is one below 200 ms falling to zero at * 500 ms. The probability of P1 is zero below 200 ms rising to * one at 500 ms and falling to zero at 800 ms. The probability * of P2 is zero below 500 ms, rising to one above 800 ms. */ up->status &= ~DGATE; if (dpulse < (200 * MS)) { p0 = 1; } else if (dpulse < 500 * MS) { dpulse -= 200 * MS; p1 = dpulse / (300 * MS); p0 = 1 - p1; } else if (dpulse < 800 * MS) { dpulse -= 500 * MS; p2 = dpulse / (300 * MS); p1 = 1 - p2; } else { p2 = 1; } /* * The ouput is a metric that ranges from -1 (P0), to +1 (P1) * scaled for convenience. An output of zero represents an * erasure, either because of a data error or pulse length * greater than 500 ms. At the moment, we don't use P2. */ return ((p1 - p0) * MAXSIG); } /* * 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[80]; /* monitor buffer */ int mldigit; /* max likelihood digit */ int diff; /* decoding difference */ int i, j; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; /* * Correlate digit vector with each BCD coefficient vector. If * any BCD digit bit is bad, consider all bits a miss. */ mldigit = 0; topmax = nxtmax = -MAXSIG; for (i = 0; tab[i][0] != 0; i++) { acc = 0; for (j = 0; j < 4; j++) { if (!(up->status & BGATE)) 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->mldigit = mldigit; vp->digprb = topmax; vp->digsnr = wwv_snr(topmax, nxtmax); /* * The maximum likelihood digit is compared with the current * clock digit. The difference represents the decoding phase * error. If the digit probability and likelihood are good and * the difference stays the same for a number of comparisons, * the clock digit is reset to the maximum likelihood digit. */ diff = mldigit - vp->digit; if (diff < 0) diff += vp->radix; if (diff != vp->phase) { vp->phase = diff; vp->count = 0; } if (vp->digprb < BTHR || vp->digsnr < BSNR) { vp->count = 0; up->alarm |= 1 << SYMERR; } else if (vp->count < BCMP) { if (!(up->status & INSYNC)) { vp->phase = 0; vp->digit = mldigit; } vp->count++; } else { vp->phase = 0; vp->digit = mldigit; up->digcnt++; } if (vp->digit != mldigit) up->alarm |= 1 << DECERR; if (!(up->status & INSYNC)) { sprintf(tbuf, "wwv4 %2d %04x %5.0f %2d %d %d %d %d %5.0f %5.1f", up->rsec, up->status, up->epomax, vp->radix, vp->digit, vp->mldigit, vp->phase, vp->count, vp->digprb, vp->digsnr); if (pp->sloppyclockflag & CLK_FLAG4) record_clock_stats(&peer->srcadr, tbuf); #ifdef DEBUG if (debug) printf("%s\n", tbuf); #endif /* DEBUG */ } up->status &= ~BGATE; } /* * wwv_tsec - transmitter second processing * * This routine is called at the end of the transmitter second. It * implements a state machine that advances the logical clock subject to * the funny rules that govern the conventional clock and calendar. Note * that carries from the least significant (minutes) digit are inhibited * until that digit is synchronized. */ static void wwv_tsec( struct wwvunit *up /* driver structure pointer */ ) { int minute, day, isleap; int temp; up->tsec++; if (up->tsec < 60 || up->status & LEPSEC) return; up->tsec = 0; /* * Advance minute unit of the day. If the minute unit is not * synchronized, go no further. */ 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]); /* * Decode the current minute and day. Set the leap second enable * bit on the last minute of 30 June and 31 December. */ 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; isleap = (up->decvec[YR].digit & 0x3) == 0; if (minute == 1439 && (day == (isleap ? 182 : 183) || day == (isleap ? 365 : 366)) && up->misc & SECWAR) 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]); /* * 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) 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 - * zero if a carry occured. 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; dp->digit++; /* advance clock digit */ if (dp->digit == dp->radix) { /* modulo radix */ dp->digit = 0; } temp = dp->like[dp->radix - 1]; /* rotate likelihood vector */ for (j = dp->radix - 1; j > 0; j--) dp->like[j] = dp->like[j - 1]; dp->like[0] = temp; return (dp->digit); } /* * 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. */ 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 * * Assuming the radio can be tuned by this program, it actually appears * as a 10-channel receiver, one channel for each of WWV and WWVH on * each of five frequencies. 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 channel, in order * to pick up minute sync and data pulses. The search for WWV and WWVH * stations operates simultaneously, with WWV on 1000 Hz and WWVH on * 1200 Hz. The probe channel rotates for each minute over the five * frequencies. At the end of each rotation, this routine mitigates over * all channels and chooses the best frequency and station. */ static void wwv_newchan( struct peer *peer /* peer structure pointer */ ) { struct refclockproc *pp; struct wwvunit *up; struct chan *cp; struct sync *sp, *rp; int rank; int i, j; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; /* * Reset the matched filter selector and station pointer to * avoid fooling around should we lose this game. */ up->sptr = 0; up->status &= ~(SELV | SELH); /* * Search all five station pairs looking for the station with * the maximum compare counter. Ties go to the highest frequency * and then to WWV. */ j = 0; sp = (struct sync *)0; rank = 0; for (i = 0; i < NCHAN; i++) { cp = &up->mitig[i]; rp = &cp->wwvh; if (rp->count >= rank) { sp = rp; rank = rp->count; j = i; } rp = &cp->wwv; if (rp->count >= rank) { sp = rp; rank = rp->count; j = i; } } /* * If we find a station, continue to track it. If not, X marks * the spot and we wait for better ions. */ if (rank > 0) { up->dchan = j; up->sptr = sp; up->status |= sp->select & (SELV | SELH); memcpy((char *)&pp->refid, sp->refid, 4); if (peer->stratum <= 1) memcpy((char *)&peer->refid, sp->refid, 4); wwv_qsy(peer, up->dchan); } } /* * 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 */ ) { struct refclockproc *pp; struct wwvunit *up; int rval = 0; pp = peer->procptr; up = (struct wwvunit *)pp->unitptr; up->mitig[up->achan].gain = up->gain; #ifdef ICOM if (up->fd_icom > 0) rval = icom_freq(up->fd_icom, peer->ttlmax & 0x7f, qsy[chan]); #endif /* ICOM */ up->achan = chan; up->gain = up->mitig[up->achan].gain; return (rval); } /* * timecode - assemble timecode string and length * * Prettytime format - similar to Spectracom * * sq yy ddd hh:mm:ss.fff ld dut lset agc stn comp errs freq avgt * * s sync indicator ('?' or ' ') * q quality character (hex 0-F) * yyyy year of century * ddd day of year * hh hour of day * mm minute of hour * ss minute of hour * fff millisecond of second * l leap second warning ' ' or 'L' * d DST state 'S', 'D', 'I', or 'O' * dut DUT sign and magnitude in deciseconds * lset minutes since last clock update * agc audio gain (0-255) * iden station identifier (station and frequency) * comp minute sync compare counter * errs bit errors in last minute * freq frequency offset (PPM) * avgt averaging time (s) */ static int timecode( struct wwvunit *up, /* driver structure pointer */ char *ptr /* target string */ ) { struct sync *sp; int year, day, hour, minute, second, frac, dut; char synchar, qual, leapchar, dst; char cptr[50]; /* * Common fixed-format fields */ synchar = (up->status & INSYNC) ? ' ' : '?'; qual = 0; if (up->alarm & (3 << DECERR)) qual |= 0x1; if (up->alarm & (3 << SYMERR)) qual |= 0x2; if (up->alarm & (3 << MODERR)) qual |= 0x4; if (up->alarm & (3 << SYNERR)) qual |= 0x8; year = up->decvec[7].digit + up->decvec[7].digit * 10; if (year < UTCYEAR) year += 2000; else year += 1900; day = up->decvec[4].digit + up->decvec[5].digit * 10 + up->decvec[6].digit * 100; hour = up->decvec[2].digit + up->decvec[3].digit * 10; minute = up->decvec[0].digit + up->decvec[1].digit * 10; second = up->tsec; frac = (up->tphase * 1000) / SECOND; leapchar = (up->misc & SECWAR) ? 'L' : ' '; dst = dstcod[(up->misc >> 4) & 0x3]; dut = up->misc & 0x7; if (!(up->misc & DUTS)) dut = -dut; sprintf(ptr, "%c%1X", synchar, qual); sprintf(cptr, " %4d %03d %02d:%02d:%02d.%.03d %c%c %+d", year, day, hour, minute, second, frac, leapchar, dst, dut); strcat(ptr, cptr); /* * Specific variable-format fields */ sp = up->sptr; if (sp != 0) sprintf(cptr, " %d %d %s %d %d %.1f %d", up->minset, up->mitig[up->dchan].gain, sp->ident, sp->count, up->errcnt, up->freq / SECOND * 1e6, MINAVG << up->avgint); else sprintf(cptr, " %d %d X 0 %d %.1f %d", up->minset, up->mitig[up->dchan].gain, up->errcnt, up->freq / SECOND * 1e6, MINAVG << up->avgint); strcat(ptr, cptr); return (strlen(ptr)); } /* * wwv_gain - adjust codec gain * * This routine is called once each second. If the signal envelope * amplitude is too low, the codec gain is bumped up by four units; if * too high, it is bumped down. The decoder is relatively insensitive to * amplitude, so this crudity works just fine. The input port is set and * the error flag is cleared, mostly to be ornery. */ static void wwv_gain( struct peer *peer /* peer structure pointer */ ) { struct refclockproc *pp; struct wwvunit *up; pp = peer->procptr; up = (struct wwvunit *)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 > 255) up->gain = 255; } else if (up->clipcnt > SECOND / 100) { up->gain -= 4; if (up->gain < 0) up->gain = 0; } audio_gain(up->gain, up->port); up->clipcnt = 0; } #else int refclock_wwv_bs; #endif /* REFCLOCK */