Precision Radio Clock for Station WWV/H Transmissions


This is a program for the Texas Instruments TMS320C25 Digital Signal Processor and Tucson Amateur Packet Radio (TAPR) DSP-93 modem. When used in conjunction with a shortwave receiver and signals transmitted by frequency and time stations WWV and WWVH., it implements a radio clock with nominal timing error less than 125 us when tracking one of the stations and frequency drift less than 0.5 PPM when not tracking either station. The clock produces an ASCII timecode that can be used to set the time of another device, such as a computer, as well as precision reference signals that can be used for other purposes, such as to drive laboratory test equipment.

The primary motivation for the program is as an example and case study of optimum receiver design using a maximum likelihood approach and matched filter, synchronous detection and soft decision principles. The clock discipline is modelled as a Markov process, with probabilistic state transitions corresponding to a conventional time-of-century clock and the probabilities of received decimal digits. The result is a performance level which results in very high accuracy and reliability, even under conditions when the one-minute beep of the signal, normally its most prominent feature, cannot be detected by ear with the shortwave receiver.

The WWV signal format is described in NIST Special Publication 432 (Revised 1990). It consists of three elements, a 5-ms, 1000-Hz pulse, which occurs at the beginning of each second, a 800-ms, 1000-Hz pulse, which occurs at the beginning of each minute, and a pulse-width modulated 100-Hz subcarrier for the data bits, one bit per second. The WWVH format is identical, except that the 1000-Hz pulses are sent at 1200 Hz. Each minute encodes nine BCD digits for the time of century plus seven bits for the daylight savings time (DST) indicator, leap warning indicator and decisecond-UTC (DUT) correction.

Program Architecture

The analog audio signal from the shortwave radio is sampled at 8000 Hz and converted to digital representation. The 1000/1200-Hz pulses and 100-Hz subcarrier are first separated using two FIR filters, a 400-Hz bandpass filter centered on 1100 Hz and a 150-Hz lowpass filter. The 1-s sync pulse is extracted using a 5-ms coherent matched filter and 8000-stage comb filter. The maximum over all stages of the comb filter establishes the first sample of the second. The 1-m sync pulse is extracted using a 800-ms noncoherent integrator and 60-stage comb filter. The maximum over all stages of the comb filter establishes the first second of the minute.

The phase of the 100-Hz subcarrier relative to the 1-s sync pulse is fixed at the transmitter; however, the audio highpass filter in most radios affects the phase response at 100 Hz in unpredictable ways. The program adjusts for each radio using two 170-ms coherent matched filters. The I (in-phase) filter is used to demodulate the subcarrier envelope, while the Q (quadrature-phase) filter is used in a tracking loop to discipline the demodulator phase.

The data bit probabilities are determined from the subcarrier envelope using a threshold-corrected slicer. The envelope amplitude 30 ms from the beginning of the second establishes the minimum (noise floor) value, while the amplitude 200 ms from the beginning establishes the maximum (signal peak) value. The slice level is midway between these two values. The negative-going envelope transition at the slice level establishes the length of the data pulse, which in turn establishes probability estimates for each of three signals, binary zero (P0), binary one (P1) and a position identifier (P2). The probability values are established by linear interpolation between the pulse lengths for P0 (300 ms), P1 (500 ms) and P2 (800 ms) so that the sum of all three values is equal to one. If the program has not synchronized to the transmitter signal, or if the data bit amplitude, signal/noise ratio (SNR) or length are less than prespecified thresholds, the data bit is considered invalid and all three probabilities are set to zero.

The difference between the P1 and P0 probabilities for each data bit is exponentially averaged in a set of 60 accumulators, one for each second, to determine the semi-static timecode data, such as DST indicator, leap second warning and DUT correction. In this design, an average value larger than a positive threshold is interpreted as one and a value smaller than a negative threshold as zero. Values between the two thresholds, which can occur due to signal fades or loss of signal, are interpreted as "don't care," and result in no change of indication.

The BCD digit in each digit position of the timecode is represented as four data bits, all of which must be valid for the digit itself to be considered valid. If so, the bits are correlated with the bits corresponding to each of the valid decimal digits in this position. If the digit is invalid, the correlated value for all digits in this position is assumed zero. In either case, the values for all digits are exponentially averaged in a likelihood vector associated with this position. The digit associated with the maximum over all of the averaged values then becomes the maximum likelihood selection for this position.

The decoding matrix contains nine row vectors, one for each timecode digit position. Each row vector includes the likelihood vector that determines the maximum likelihood decimal digit for that digit position, along with other related data. The maximum likelihood digit for each of the nine digit positions becomes the maximum likelihood time of the century. A built-in transition function implements a conventional clock with decimal digits that count the minutes, hours, days and years, as corrected for leap seconds and leap years. The counting operation also rotates the likelihood vector corresponding to each decimal digit as it advances. Thus, once the seconds and minutes have been synchronized, the maximum likelihood time in any minute should correspond to the BCD timecode transmitted in that minute.

Each row of the decoding matrix includes the likelihood vector, clock digit, a compare counter and the difference (mod 10) between the current clock digit and most recently determined maximum likelihood digit. If a digit probability exceeds the decision level and the difference is constant for a number of successive minutes in any row, the maximum likelihood digit replaces the clock digit in that row. When the differences are zero for all rows, the clock is synchronized and delivers correct time to the integral second. The millisecond within the second is derived from the A/D sample clock, which runs at 8000 Hz and drives all system timing functions.

The A/D sample clock is derived from the DSP-93 CPU clock and a system of hardware counters. Its frequency is disciplined by a frequency-lock loop (FLL) which operates independently of the data recovery functions. At averaging intervals determined by the measured jitter, the frequency error is calculated as the difference between the most recent and the current epoch divided by the interval. The sample clock frequency is then corrected by this amount using an exponential average. When first started, the frequency averaging interval is eight seconds, in order to compensate for DSP-93 CPU clock frequency offsets up to 62 PPM. Under most conditions, the averaging interval doubles in stages from eight seconds to over 1000 seconds, which results in an ultimate frequency precision of 0.125 PPM, or about 11 ms/day. This is consistent with the stability of a TCXO clock oscillator and somewhat better than the uncompensated oscillator used in the DSP-93.

Measured Performance

It is the intent of the design that the accuracy and stability of the indicated time be limited only by the characteristics of the propagation medium. Conventional wisdom is that synchronization via the HF medium is good only to a millisecond under the best propagation conditions. The performance of the program is clearly better than this, even under marginal conditions. Ordinarily, with marginal to good signals and a frequency averaging interval of 1024 s, the frequency is stabilized within 0.1 PPM and the time within 125 us. The frequency stability characteristic is highly important, since the clock may have to free-run for several hours before reacquiring the WWV/H signal.

The expected accuracy over a typical day has been determined using an oscilloscope and cesium oscillator calibrated with a GPS timing receiver. With marginal signals and allowing 15 minutes for initial synchronization and frequency compensation, the time accuracy determined from the WWV/H 1-s sync pulse is reliably within 125 us. In the particular DSP-93 used for program development, the uncorrected CPU clock frequency offset is 45.8+-0.1 PPM. Over the first hour after initial synchronization, the clock frequency drifts about 1 PPM as the frequency averaging interval increases to the maximum 1024 s. Once reaching the maximum, the frequency wanders over the day up to 1 PPM, but it is not clear whether this is due to the stability of the DSP-93 clock oscillator or the changing height of the ionosphere. Once the frequency has stabilized and after loss of the WWV/H signal, the frequency drift is less than 0.5 PPM, which is equivalent to 1.8 ms/h or 43 ms/d. This can result in a step phase correction up to several milliseconds when the signal returns.

The measured propagation delay from the WWV transmitter at Boulder, CO, to the receiver at Newark, DE, is 23.5+-0.1 ms. This is measured to the peak of the pulse after the 1-s sync comb filter and includes components due to the ionospheric propagation delay, communications receiver delay and DSP-93 delay. The propagation delay can be expected to change about 0.2 ms over the day, as the result of changing F-layer height. The measured DSP-93 delay is 5.5 ms, most of which is due to the 400-Hz bandpass filter and 5-ms matched filter.

Program Operation

The program begins operation immediately upon startup. It first searches for 1-s sync, then 1-m sync, which can take up to several minutes, depending on signal quality. When 1-s sync has been acquired, LED4 on the DSP-93 front panel begins to flash. When 1-m sync has been acquired, LED3 begins to flash as well, with the duration of the flash following the pulse-width modulated data pulse. At this time, the tracking loop begins to discipline the demodulator phase to nominal zero offset relative to the subcarrier phase, which may take several minutes. LED5 is on when the data pulse has sufficient amplitude, SNR and duration; however, under marginal conditions or loss of signal, LED5 may be on only occasionally throughout the minute.

The program then begins accumulating likelihood values for each of the nine digits of the clock, plus the seven miscellaneous bits included in the WWV/H transmission format. LED6 flashes when the likelihood value is below the decision threshold for a digit or miscellaneous bit. LED7 flashes when the current clock digit disagrees with the maximum likelihood digit. The minute units digit is decoded first and, when several repetitions have compared correctly, the remaining eight digits are decoded. When several repetitions of all nine digits have decoded correctly, which normally takes 15 minutes with good signals and up to an hour when buried in noise, the clock is set within 125 us and can be read within 500 μs with the ASCII timecode formats described later.

At the same time, the clock discipline algorithm refines the frequency offset correction for use during times when the WWV/H signal is unavailable. The algorithm refines the offset using increasingly longer averaging intervals to 1024 s, where the precision is about 0.1 PPM. With good signals, it takes well over two hours to reach this degree of precision; however, it can take many more hours than this in case of marginal signals.

To work well, this program needs a communications receiver with good audio response at 100 Hz. Most shortwave and communications receivers roll off the audio response below 250 Hz, so this can be a problem, especially with receivers using DSP technology, since DSP filters can have very fast rolloff outside the passband. Some DSP transceivers, in particular the ICOM 775, have a programmable low frequency cutoff which can be set as low as 80 Hz. However, this particular radio has a strong low frequency buzz at about 10 Hz which appears in the audio output and can affect data recovery under marginal conditions. Although not tested, it would seem very likely that a cheap shortwave receiver could function just as well as an expensive communications receiver.

Front Panel Indicators

The status of the Program is shown by the DSP-93 front panel LEDs, which are interpreted as follows:
Signal level indicator. Adjust the input signal level until this LED flashes on occasional peaks.
Activity indicator. Intensity corresponds to the fraction of idle processor time; i.e., not doing signal processing functions.
Data indicator. Follows the 100-Hz subcarrier modulation that carries the data bits. This is active only when in 1-s sync.
Sync indicator. Blinks to indicate 1-s and 1-m sync. Short blinks follow the 1-s sync pulse, while long blinks follow the 1-m sync pulse. This is active in 1-s sync before the clock has synchronized when and always after the clock has synchronized. If the clock has not been in both 1-s sync and 1-m sync during any minute, status bit 0 is set in the following minute.
Valid data. Turned on for one second only when in 1-s and 1-m sync and following a valid data bit. An invalid data bit has assumed bit probability of zero, but does not necessarily indicate a clock error. Note that there is no data bit in second zero of the minute. If there too many invalid data bits during any minute, status bit 1 is set in the following minute.
Low decoded signal probability. On only in 1-s sync and when a correlated bit or digit probability is below the decision threshold. The associated bit probability is zero, but does not necessarily indicate a clock error. LED6 when on causes status bit 2 to be set in the following minute.
Decoding error. On only in 1-s sync and when the maximum likelihood digit differs from the current clock digit. The associated digit is not used to set the clock, but does not necessarily indicate a clock error. LED7 when on causes status bit 3 to be set in the following minute.
DSP-93 Power indicator

Program Commands

The DSP-93 program responds to a number of commands, which consist of a single case-insensitive letter followed in some cases by one or more digits. Control characters, including CR and LF, are ignored, except in the C and H commands. The station to be monitored is selected using the C (WWV) or H (WWVH) commands. The only affect these commands have is to select which 5-ms matched filter to use for the 1-s sync, which tone (1000 Hz or 1200 Hz) to use for the 1-m sync integrator, and which propagation delay value to use. These commands can be followed by a decimal number, which specifies the propagation delay in 125-us units, and then a CR to enter the value. A C or H command followed immediately by a CR selects the station, but does not change the propagation delay.
This command selects the UART baud rate in bps. The digit argument is interpreted as follows:

0            300
1     600
2     1200
3     2400
4     4800
5     9600
6     19200
7     38400
8     76800
9     153600
The default is set by the program loader, which is 19200 with the TAPR firmware monitor.
This command enables or disables either the RTS or DTR line from the UART to follow LED4. Note that internally the UART RTS output is connected to the CTS pin (8) on the DB9 connector. The digit argument is interpreted as follows:

0     RTS and DTR always high
1     RTS follows LED4, DTR always high
2     DTR follows LED4, RTS always high.
This command resets the program to the initial state at the point control is transferred by the TAPR firmware monitor to the program. S{0-9} This command sets the nine digits of the clock from low to high order. This command is useful in order to verify proper rollover of the leap second and leap year, for example.
This command returns the timecode string in the following format preceded by a CR/LF sequence. This format conforms to Format 2 as used by the Spectracom WWVB and GPS radio clocks and supported by the current version of the Network Time Protocol (NTP) software for Unix, Windows and VMS.

sqyy ddd hh:mm:ss.fff ld
s     sync indicator (? or SP)
q     quality character (SP)
yy    year of century
ddd   day of year
hh    hour of day
mm    minute of hour
ss    second of minute
fff   millisecond of second
l     leap second warning SP or L
d     DST state S, D, I, or O

The on-time reference is the start bit of the first (CR) character. The sync indicator s is '?' until the clock is set by valid received data and SP after that. The quality character q is always SP in this format.

The leap second warning l is L if a leap second is expected at the end of the month and SP if not. By international agreement, the leap second is inserted in the last second of the last day of June or December and is numbered second 60 of that minute. The WWV/H timecode progression just before, during and after the leap on 30 June 1997 is shown below, where the first column is the year, followed by the day, time, leap second warning, DST state and DUT correction (see U command):

97 181 23:59:59 LD -5
97 181 23:59:60 LD -5
97 182 00:00:00 D  +5

Note the seconds progression includes second 60 of the last minute and the day number advances at the next second. Also, the leap warning disappears and the DUT correction -5 changes to +5 at that time; however, in actual use these changes are delayed a few minutes as the result of the integration functions.

The DST state is determined from the DST1 and DST2 bits of the WWV/H timecode. The DST1 bit is set to one at 0h UTC on the day of transition to daylight time and thereafter. The bit is set to zero at 0h UTC on the day of transition to standard time and thereafter. The DST2 bit bit follows the state of the DST1 bit, but delayed 24 hours later. The values are decoded as follows: S standard time D daylight time I between 0h on the day to daylight time until 0h the next day O between 0h on the day to standard time until 0h the next day.

This command returns the timecode string in the following format followed by a CR/LF sequence. This format is intended for experimentation and evaluation and may be changed in future.

sq a yy ddd hh:mm:ss.fff ld dut bbbbb fffff ggggg uuuuu

where the elements of the timecode are as described for the T command, with the addition of the following:

q     quality character
a     station select C (WWV) or H (WWVH)
dut   DUT sign and magnitude in deciseconds (tenths of a second)
bbbbb count of invalid data bits during the current minute
fffff frequency offset in units of about 1.9 ns/s
ggggg frequency averaging interval in seconds
uuuuu interval since last clock update in minutes

The on-time reference is the start bit of the first (s) character. The quality character q encodes four bits of status formatted in hex. Reading from high to low order, the bits are:

0     low 1-s or 1-m sync pulse amplitude
1     too many invalid data bit errors
2     low decoded digit probability
3     timecode digit decoding error

The DUT sign and magnitude encode the UT1 time correction in deciseconds. The UT1 correction is used by astronomers and navigators to establish a precise star transit time.

Debugging and Monitoring Commands

The program includes a primitive toolkit useful for debugging and evaluation. Debugging mode can be turned on with one of the D commands and off with the D0 command. There are four debug formats selected by the D1, D2, D3 and D4 commands, respectively. The output is inclusive; that is, one, two or all three formats can be enabled at the same time, but all are disabled by the D0 command.
This command produces one line at second 59 of each minute in the same format as the U command described above.
This command produces one line for each decoded timecode digit in thefollowing format:

c m d n q sssss aaaaa mmmmm ddddd nnnnn eeeee qqqqq

c     clock digit
m     maximum likelihood digit
d     difference (mod 10) between clock and maximum likelihood digits
n     compare counter (max 3)
q     quality character
sssss 1-s sync epoch
aaaaa 1-s sync amplitude (threshold 1500)
mmmmm 1-m sync amplitude (threshold 500)
ddddd 100-Hz subcarrier amplitude (threshold 1000)
nnnnn 100-Hz subcarrier SNR (dB times 10 - threshold 100)
eeeee digit correlator amplitude (threshold 1000)
qqqqq digit correlator SNR (dB times 10 - threshold 30)

This command produces one line each second in the following format:

aaaaa eeeee ppppp fffff ccccc jjjjj vvvvv

aaaaa 1-s sync amplitude
eeeee median filter epoch of 1-s sync
ppppp median filter time offset
fffff time offset since last frequency update
ccccc compare counter (max 10)
jjjjj jitter counter (max equal to averaging cycle)
vvvvv averaging cycle counter

This command produces one line each second in the following format:

aaaaa eeeee ppppp fffff ccccc

aaaaa Second number in the minute
eeeee data bit envelope amplitude
ppppp data bit envelope SNR (dB times 10)
fffff raw data bit pulse length (relative)
ccccc integrated data bit pulse length (relative)

This command produces one line at second 58 of each minute in the following format:

ddddd eeeee fffff ggggg hhhhh iiiii jjjjj

ddddd DST1 integrator value
eeeee DST2 integrator value
fffff leap warning integrator value
ggggg DUT+- integrator value
hhhhh DUT1 integrator value
iiiii DUT2 integrator value
jjjjj DUT4 integrator value

Signal Generation and Monitoring

This command selects certain signals to be monitoring in real time with an oscilloscope connected to the D/A converter. The synthetic signals are synchronized with respect to the WWV/WWVH transmitter, so may be useful as precision time and frequency sources for other devices, such as laboratory test equipment. Note that the stuffing cycles necessary to produce the standard frequencies from the DSP-93 hardware sample clock produce some phase modulation in the signals, which appear as a low frequency buzz on the 1000/1200-Hz signal. The signals are selected with the O{1-9} command, where the digit is:

1    input signal (after peak clipping, but before other processing)
2    output of 400-Hz bandpass filter
3    output of 150-Hz lowpass filter
4    output of 1-s sync comb filter
5    output of matched filter I channel
6    output of matched filter Q channel
7    synthesized one-second ramp in 125-us increments
8    synthesized 1000-Hz (WWV) or 1200-Hz (WWVH) sinewave
9    synthesized 100-Hz sinewave
0    data bit probability - positive (1), negative (0)
For the DC-coupled signals, the oscilloscope should be connected to JP210, pin 2; that is, the pin closest to the power regulators. In the case of the ramps and bit probability, the lowpass filter in the DSP-93 D/A converter limits the risetime and introduces a moderate degree of ringing.

A speaker can be connected to the DSP-93 output to monitor the operation of the program. The O1 signal is useful for general watchkeeping. The O2 signal is probably the best indication that the WWV/H signal is present. If a faint tick can be heard, the program can probably synchronize to the signal. Normally, at the decision threshold where the tick can barely be heard, the WWV/H signal cannot be heard in the unprocessed receiver output, although a faint beep may be heard marking the first second of the minute. The O3 signal is a good indication of the subcarrier SNR; however, there can be a huge degree of integration in this signal and it may take many minutes to stabilize. The O4 signal is best viewed with an oscilloscope, where the waveform should appear as a trapezoid rising at the second and decaying either 200, 500 or 800 ms later. When viewed on an oscilloscope with a calibrated sweep delay and triggered by a PPS signal from a source calibrated to UTC, the O5 signal becomes a precision proof-of-performance indicator.

Operation Notes

The particular algorithms implemented in the program have been thoroughly tested over many months and seasons, but the various thresholds have not been experimentally verified under all conditions. This applies especially to the decision thresholds, called bottom fishers, which determine whether a particular signal is usable or not. Since the program is designed to be very aggressive in sifting signals from noise, if these thresholds are set too low, the program may decode noise as UTC. On the other hand, if these thresholds are set too high, usable signals may be ignored and the clock drift away from UTC. The current thresholds are defined in the output formats described above and in the source code. They are probably set too low in an effort to scrape the weakest signals from the noise. In normal operation, the values are several times the thresholds.

There is a necessary design compromise between weak-signal sensitivity and DSP-93 CPU oscillator frequency tolerance. When first started, the program must discipline the sample clock to within 12.5 PPM, in order to reliably identify the 1-s sync signal epoch. At the smallest averaging interval (8 s) and largest jitter threshold (500 us) which can reliably extract the epoch, the tolerance must be less than 62.5 PPM. As the frequency offset of the DSP-93 used for program development is over 45 PPM, this could be a problem in some units. While not tested, the program should cope with tolerances greater than 62.5 PPM, but it will take significantly longer to synchronize.

Following is a list of shortcomings and proposed enhancements, which may be fixed/implemented in future:

  1. The program should provide optional local timezone correction, standard/daylight time adjustment and DUT correction.
  2. The program should provide optional formats emulating other radio clocks on the market.
  3. The program should provide port selection, gain control and possibly other features useful for system integration.
  4. The program should provide automatic UART baud rate control using the AT command.
  5. The program should be rewritten in portable C for use on a fast workstation with integrated audio codec.
  6. The program should support the Canadian time/frequency station CHU. This would require provision for a 300-bps modem and not much else.


Written by David L. Mills, W3HCF; this update 26 May 1998


This program is a work in progress. Your comments and suggestions are most welcome.

David L. Mills (