[time-nuts] DIY Loran-C frequency receiver update (2)

[Poul-Henning Kamp]

In message <12050.1222638372 at critter.freebsd.dk>, Poul-Henning Kamp writes:

>Some time ago I send an email about my little project to create a
>LORAN-C frequency receiver for time-nuts purposes, here is an update.

I have gotten quite far with the software now, it can pull in any
signal that can reasonably be expected to be heard at my location.

The software tries to latch onto the strongest loran-station in the
GRI, but due to the approximately equal distance to the stations
it waffles between three different stations in the 9007 chain.

Here are the six signals it can dredge out of my test data:

	http://phk.freebsd.dk/misc/loran_20091013.pdf

The stations are:

	Signal	Location	Relativ vol.	Time to lock
	7499M	Sylt		1.000		16 s
	6731Z	Sylt		0.732		15 s
	9007M	Eiðe		0.043		98 s
	9007X	Bø		0.036		69 s
	7001M	Bø		0.022		60 s
	9007Y	Værlandet	0.021		90 s

If the software was told to skip the stong station and take one
of the weaker, more stations could be tracked, just like increasing
the integration time would enable pulling weaker stations out
of the noise.

Neither of this is the point however, the idea is to measure
the frequency of the strongest signal available.


The way the software works is the following:

Phase 1:

Pick out every 100th sample (one per msec) and average these over
the GRI until a pattern of four spikes stand out in either ^v^^
(master) or ^^^v (slave) pattern, at least 2.5 standard deviations
over the set of up to 1000 buckets.

This tells us where the signal we want is, and which code it uses.

The first two plots from the top shows the results of this, with
the red bars indicating the best candidate signal.

Phase 2:

Integrate the first 0..500 and 1000..1500 samples of each gri
and multiply the first set with the second set.  Once the
relative result goes above or below +/- 0.75 we know if we
have the A or B period of the signal.  This is the same for
master and slave codes, since they both have "++..." and "+-..."
in their A and B periods.

The middle plot shows the 1000 buckets and the result of this.

Phase 3:

Integrate the 0..1000 samples summed over the full code.

This is where we will identify the 3rd zerocrossing and track
its position relative to the sample clock.

The fourth plot shows this.

The fifth plot is a debugging aid and shows the full FRI
accumulated along with the fourth plot, at the left is the
A code of the tracked signal, in the middle the B code and
for some of the signals, other members of the chain are visible
along the line.  See for instance page 2 where all the 6731
stations are visible.

Phase 4: is not coded yet, and consists of identifying the
3rd zero crossing and interpolating it's position between
samples.


One intesting detail is that the software does not work for
the 8000 Chayka chain, because averaging over the FRI or GRI
does not filter out broadcasting stations transmitting on
whole multiples of kHz.

The 8000 sample average corresponds to 8 msec of spectrum, so a
signal of 198 kHz would fit exactly 198 * 8 = 1584 periods into
the window and these will integrate just nicely.

Interestingly, this also works for pulling out frequency: integrating
the 8000 bucket multiplied by sin and cos for that frequency gives
a nice phase-vector for the BBC longwave station, or any other
station running on a multiple of 1 kHz.

This should not be a problem for north american chains which
have the XXX0 GRI format, but the averaging algorithm may
not be as effective on these as on the XXXX european GRIs.

Ohh, and the code is 592 lines of C code (exclusive the functions
to read data from my compressed files and to emit postscript plots).

-- 
Poul-Henning Kamp       | UNIX since Zilog Zeus 3.20
phk at FreeBSD.ORG         | TCP/IP since RFC 956
FreeBSD committer       | BSD since 4.3-tahoe    
Never attribute to malice what can adequately be explained by incompetence.