Wednesday, October 31, 2018

Decoding STANAG 4285 using KiwiSDR and GNURadio

Below are some examples of STANAG 4285 signals on shortwave picked up by a KiwiSDR located in Scandinavia. These signals were decoded live using the new KiwiSDR GNURadio block gr-kiwisdr, see this blog post, and another GNURadio block I started to work on, called gr-digitalhf:
  • gr-digitalhf contains an adaptive decision feedback equalizer (DFE) and a local oscillator for Doppler correction, implemented in C++. These are separated from the description of the physical layer (synch preamble, frame format, constellation, ...) which is provided in Python. You may have a look at the description of STANAG 4285 here which is an almost literal translation of the standard document.
  • The DFE unscrambles a given symbol, computes a hard decision using the constellation valid for the unscrambled symbol, and then uses the scrambled symbol from the hard decision for computing the error. For this it only needs to have (1) an estimate of the Doppler (=frequency offset) correction and it needs to know (2) if a given symbol is known or unknown. Both, (1) and (2) are provided by the description of the physical layer in Python.
  • Currently the implementation is limited to providing soft decisions. The two final steps, deinterleaving and convolutional decoding, are still missing.
The aim of gr-digitalhf is to make it easy to adapt it to other digital modes based on N-ary PSK/QAM modulation which have a more complicated physical layer compared to STANAG 4285, like, e.g., MIL-STD-188-110A.

STANAG 4285 signal (BPSK unscrambled data, 600 baud) with severe multi-path.
Another STANAG 4285 signal (BPSK unscrambled data, 600 baud)
An example of a 1200 baud STANAG 4285 signal (QPSK unscrambled data)


Monday, October 22, 2018

KiwiSDR GNURadio block

A first version of a GNURadio interface gr-kiwisdr for the KiwiSDR is available on GitHub. It consists of a GNURadio source block which connects to a given KiwiSDR and provides a stream of IQ samples.

At this moment gr-kiwisdr is not yet complete but fully working. Missing items include support for
  • other modes than IQ mode (not sure if this is really useful)
  • user login and non-default AGC settings
  • stream tags containing GNSS timestamps

For connecting to the KiwiSDR websocket server boost.beast is used. Because the available boost versions, e.g., on Ubuntu, are fairly old and boost.beast is a recent addition to boost, gr-kiwisdr contains git submodules for boost.asio and boost.beast (both are header-only libraries).

Any feed-back is very much welcome.
As an example, the preamble cross-correlation of a STANAG 4285 signal is shown below. The two secondary cross-correlation peaks are there because the 80-symbol long STANAG 4285 preamble consists of repetitions of a 31-bit long pseudo-random sequence generated by a LFSR.

Real-time preamble cross-correlation for a STANAG4285 signal.


GRC flow-graph for the plot shown above.

Burst signal monitoring using kiwirecorder.py

Recently burst signals on 7001.8 kHz (center) were observed on KiwiSDRs in Scandinavia. Simimlar bursts can be observed on 4089.8 kHz.

The bursts repeat every 5 seconds and are PSK8 modulated with 2400 symbols/sec. All bursts seem to be identical and have a weak peak in the auto-correlation function at ±215 samples corresponding to 43 symbols.

The bursts were recorded using kiwirecorder.py with the option --squelch-threshold 10 which starts a new file for each burst. The --squelch-tail option is new and can be used to set the time the squelch remains open after the signal is below threshold.

Tuesday, October 16, 2018

STANAG 4285 signal on 442 kHz

A STANAG 4285 signal was detected on 472 kHz on the Kongsfjord KiwiSDR on 472 kHz:

First I thought it was a fake signal but it was also present on the SM2BYC KiwiSDR.

The bitstream obtained after convolutional decoding and deinterleaving consists of frames delimited by the LFSR sequence F1, using the notation of this blog post:

TDoA maps point to somewhere on the coast of Norway:

Friday, September 28, 2018

KiwiSDR GNSS position solutions using an extended Kalman filter

This is an update to this blog post, see also https://github.com/jks-prv/Beagle_SDR_GPS/issues/147. References for using Kalman filters for GNSS positioning are, e.g., DST-Group-TR-3260.pdf and these proceedings from ION GPSGNSS 2003.

The plots below show time series of variations in ECEF X,Y,Z coordinates, the oscillator correction, and the number of available satellites.

SPP - single point position; EFK - extended Kalman filter


The process noise covariances are set to 5e-5 for the X,Y,Z coordinates (KiwiSDRs are not expected to move) and the covariances for the time and the oscillator correction are build from Allan variance parameters as described in a497248. It would be interesting to measure the Allan variance of the used oscillator.

The advantages of using an extended Kalman filter include:
  • Once a position solution has been obtained, the Kalman filter can be updated with fewer than four satellites: this will improve the availability of those KiwiSDRs with poor GNSS reception for TDoA work.
  • The residuals per link in the Kalman filter update are used to exclude satellites with bad time measurements.
  • Large deviations from the nominal position due to bad satellite geometries are avoided.

Tuesday, September 25, 2018

New version of TDoA code available

A version of the 1st idea mentioned in the last blog post is now implemented in the latest version of the TDoA code and thanks to John it is available as an option in the KiwiSDR TDoA service.

The changes in the new TDoA algorithm are the following:

(1) Sanity checks on the recorded IQ wav files


Instead of aborting processing when a single file is found to be bad, stations are now excluded. If more than 1 good station is left, the processing proceeds without the excluded stations. The following checks are performed:
  • check for recent GNSS position fixes
  • validity check on then GNSS timestamps
  • check if |IQ|!=0
  • check if the overlap of all wav files is long enough (>10 sec)

(2) Finding clusters in the cross-correlation peaks


For a given pair of stations, peaks in the cross-correlations are found as before. Then these peaks are assembles into up to four clusters using an adaptive clusterization algorithm. It starts by assuming that there are four clusters. If the four clusters are not valid, the algorithm starts again, this time assuming that there are three clusters. When all clusters are valid the algorithm terminates.

The used clusterization algorithm is iterative:
  • It starts by assuming that there are N clusters and puts N+1 boundaries bi at equidistant places between the minimum and maximum time of the cross correlation peaks.
  • Then it recomputes the boundaries as follows: Let ti,i+1 be the average of the mean of the cross correlation peaks between bi and bi+1. The boundaries are iteratively updated as bi=(ti-1,i+ti,i+1)/2 until they become constant.
  • A cluster is called valid if the distribution of entries in it is not compatible with a random distribution.
  • If all clusters are valid or if there is only one cluster (N=1) then the clusterization algorithm terminates.

(3) Selection of those clusters which optimally fulfill the consistency equations

  • If there are three stations, cross-correlation times of the same signal fulfill (12)-(13)-(23)=0.
  • When there are four stations, there are more of these consistency equations, e.g., (12)-(14)-(24)=0, (23)-(24)-(34)=0 etc.
  • For all combinations of clusters all consistency equation are evaluated. Each equation provides the deviation from 0 in terms of number of sigmas (nsigma) using the mean and RMS of the clusters.
  • In the end the combination of clusters is selected which for all consistency equations results in the smallest nsigma
  • If the nsigma value in any given equation is greater than 3, stations are iteratively excluded until 3 stations are left. The combination of excluded stations which brings max(nsigma) below 3 is not used later on for making the TDoA maps.

Thursday, August 23, 2018

New ideas for TDoA multi-lateration

(1) Finding the common propagation mode


Let's consider the case of three receivers receiving the same signal.
  • Each receiver receives the signal at times t1, t2, t3.
  • The cross-correlations produce time differences Δt1,2, Δt1,3, Δt2,3.
  • Trivially, Δt1,2 = Δt1,3 - Δt2,3 (and cyclically) (*)
This situation is illustrated in the plot below.

Common propagation mode


As a second step let's consider the situation when one of the three receivers receives the common signal via two propagation paths, and the cross-correlations to the signals received at the two other receivers have approximately the same magnitude.

Currently the TDoA software performs for each cross-correlation a clustering algorithm, and in the situation shown in the plot below it may well be that the following time differences are obtained: Δt1,2 = -50, Δt1,3 = -100, Δt2,3 = -140. Note that these do not fulfill the identity (*).

Two propagation modes for one receiver

As a consequence, a trivial stage of pre-processing is given by selecting those cross-correlation peaks which are consistent according to (*).

Update 9/22/2018


(2) Using constraints (number of degrees of freedom)


Now let us count the number of degrees of freedom for the case of N receivers receiving the same signal: there are
  • 2 unknowns, (lat,lon), for the most likely position
  • N unknown virtual heights of which only (N-1) are observable in the time differences (adding a common time offset does not change Δti,j)
  • N(N-1)/2 hyperbolæ (=constraints)
Note that instead of virtual heights we might as well add time offsets δt to the time differences, Δti,j+δtj, j≠1, and in this way become agnostic w.r.t. the specifics of ionospheric propagation.

As a summary there are N+1 unknowns and N(N-1)/2 constraints:

#receivers #unknowns #constraints #constraints-#unknowns
3 2(pos)+2(δtj, j≠1) 3 hyperbolæ -1
4 2(pos)+3(δtj, j≠1) 6 hyperbolæ 1
N 2(pos)+(N-1)(δtj, j≠1) N(N-1)/2 hyperbolæ N(N-1)/2-N-1

So for 3 receivers there is a 1-dimensional family of solutions while for N≥4 receivers the system of equations is overdetermined.  This sounds good enough to make me suspicious if it is indeed true.



It will be interesting to try these ideas on real data obtained from KiwiSDRs. If anyone spots an error in the logic above please do let me know.





When using N stations there are N-1 independent measurements contained in the time differences δtij, j>i. Two of these are needed to determine the likely position, so N-3 out of the N-1 δti, i≠1, can be found:

#receivers #unknowns #measurements #fixed/total δtj
3 2(pos)+2(δtj, j≠1) 2 time differences 0/2
4 2(pos)+3(δtj, j≠1) 3 time differences 1/3
5 2(pos)+4(δtj, j≠1) 4 time differences 2/4
6 2(pos)+5(δtj, j≠1) 5 time differences 3/5
N 2(pos)+(N-1)(δtj, j≠1) (N-1) time differences (N-3)/(N-1)

Sunday, August 5, 2018

Detection of k=7 r=1/2 convolutional codes

We start with a random binary sequence xi which is encoded using a k=7 r=1/2 convolutional decoder into
  yI = [u0,v0,u1,v1,...]  ,
with
  ui = ajxi+j and vi = bjxi+j  ,
where "+" denotes XOR and repeated indices are summed over from 1 to 7.

The property used to find a,b is the following:
  bjuj+i = ajvj+i  .   (*)

For this we start with the full state space of 127×127 values of a,b. For each pair (ui, vi) the property (*) is evaluated and states which fail it are pruned from the search space.

The screen shot below shows the output of a simple octave script demonstrating this approach where further details of this algorithm can be found. For an encoded random bit sequence, after consuming about 28 encoded bits yI, (i.e., 14 iterations), the correct a,b values are found.

Output of signal-analysis/m/test_conv.m


This method must be already known and I would be grateful if someone showed me a reference to it.

Sunday, July 15, 2018

TDoA SuperDARN/Finland

Today signals from the SuperDARN radar in Finland were active around 12375 kHz and the most likely position matches quite well with the known CUTLASS radar position near Hankasalmi.

TDoA maps


TDoA cross correlations


There are pulse trains spaced by integer multiples of 2400 μsec (0,9,12,20,22,26,27), as descibed here:

SuperDARN pulse train

30 of these pulse trains are repeated, approximately each 95.3 msec, then there is a gap of two frames and some delay, and then the next 30 pulse trains follow:

SuperDARN pulse frames

The transmissions seem to repeat each minute with a small gap at the end of the minute.

Another way to verify the used pulse sequence is to compute the autocorrelation of each line of the plot above. All differences between all combinations of the individual pulses show up as peaks in the autocorrelation:
frame number vs. autocorrelation lag

Friday, July 13, 2018

KiwiSDR TDoA (LORAN)

Below is an example of using the new KiwiSDR TDoA extension on LORAN signals, see The bias on the most likely position is quite small. Note that despite the fact than one sample at 12kHz is 83.3μs long, the RMS values of the cross correlations are ≤ 6μs. This is because the peaks in the cross correlations are fitted with a 2nd order polynomial.

Fig 1

Fig 2


Fig 3


Fig 4


Fig 5


Fig 6


Fig 7

Wednesday, July 11, 2018

Propagation delay maps based on IRI2016 electron density profiles

In order to understand better the propagation delay maps used for this blog post, the animation below was created. References for the used method are arXiv:1104.2248 and this well-known book.

Up to 2-hop propagation is considered. In case there is more than one mode of propagation, the smallest non-ground-wave propagation delay is shown in the maps. If there is no ionospheric reflection, the ground-wave delay is shown.
  • hop=0: no ionospheric reflection -> ground-wave propagation
  • hop=1: one reflection
  • hop=2: two reflections 
Propagation delays based on IRI2016 electron density profiles

Propagation delays based on IRI2016 electron density profiles

Saturday, July 7, 2018

KiwiSDR TDoA extension

KiwiSDR v1.196 includes a KiwiSDR extension (based on https://github.com/hcab14/TDoA)

Note that ionospheric delays are not taken into account, so depending on the geometry of receivers and the target the TDoA maps will be biased.

Many thanks to John Seamons, ZL/KF6VO.

In parallel I am working on a real-time version of TDoA using python/matplotlib/cartopy.

Thursday, June 21, 2018

Interesting signal (TDMA?)

Today a signal on 5099 kHz USB (5100.8 kHz center) was picked up on several KiwiSDRs in the Northwest of Europe.

This signal consists of a sequences of sharp pulses and of tones:



Zoom around the tone:



Aligning the signal in frequency to the beginning of the tone, it becomes visible that there are two tones in three intervals. It might be a coincidence but the length of the intervals is close to 6/128 seconds which may or may not be relate to the information (here or here):

It is interesting that some residual modulation is visible, i.e., the phase is not perfectly constant or linearly rising.

The next three plots show abs(z) assuming different frame lengths:
  1. 854 samples (9.11/128 sec):



  2. 844 samples (9/128 sec):



  3. 281 samples (3/128 sec):



Monday, June 11, 2018

TDoA with propagation delays based on IRI2016 electron densities

This is a follow-up of this blog post where the method referred to there is applied to TDoA maps. The first example is a STANAG 4285 signal (FUO Toulon) centered on 8436.4 kHz (USB: 8432.6 kHz)

TDoA maps using ground-wave propagation delays


TDoA using propagation delays based on IRI2016


The maps using propagation delays based on IRI2016 electron densities match better the (assumed) true position of the transmitter. It is a small effect but noticeable, and can be seen better in the side-by-side comparison below:

TDoA comparison



Another example are the TDoA maps from the last blog post where the TDoA maps based on more realistic propagation delays are more consistent between the two measurements; both indicate a position on the South coast of Cornwall:

TDoA comparison


TDoA comparison


Sunday, June 10, 2018

TDoA code update

The TDoA code now does not explicitly use gnuplot any more. All plotting is done in octave. This will it make easier to extend the code (work in progress).

Recently some interesting signals have been independently observed by several people, see, e.g., tweet. These are frequency hopping signals.

Tuning a KiwiSDR to 7812 kHz two of the signal fit into its passband. Below are TDoA maps for these signals. Note that these maps are likely biased by the fact that some of the KiwiSDRs might receive the ground-wave signal and others a sky-wave signal.

TDoA maps


TDoA cross correlations


TDoA maps


TDoA cross corrlelations

Saturday, June 9, 2018

Towards ionospheric ray-tracing

For the work on HF TDoA a realistic method for the signal delay is needed. Up to now, the great-circle distance divided by the speed of light is used which of course is incredibly naive. Nevertheless it seems to work, at least to some degree.

Rather than attempting a full-fledged ray-tracing simulation one can use the theorems of Breit and Tuve and Maryn's theorem for obtaining the signal propagation delay from virtual vertical height profiles, see, e.g, arXix:1104.2248 and here.

In order to evaluate this technique let us use the oblique ionograms measured by the University of Twente web SDR (thanks to PA3FWM for making these available! and for him and G3PLX for explanations):

The plots below show the time delay vs. UTC for two days in June at two different frequencies. Overlaid are time delays obtained by
  • VOACAP: all propagation modes with probability>0.3 are shown,
  • simulations explained above using electron density profiles obtained from the international reference ionosphere model IRI2016
(For this comparison the oblique ionograms have been corrected for a time offset of ~0.3ms whose origin is not yet clear)

All plots are generated using octave and this octave binding of the IRI2016 model.

As the IRI model is a climatological model (as is VOACAP), the agreement between model calculation and data is not expected to be perfect. Only the general features of these backscatter ionograms are reproduced, such as the 2-hop E-layer reflection for lower frequencies during mid-day but, e.g., not sporatic E-layer reflections in the late afternoon/evenings.

Note also that this simple method does not take into account magneto-ionic effects which become important for lower frequencies (and at local night-times).







Wednesday, June 6, 2018

Blind detection of checksums (1)

This is a generalization of T207_test.m to arbitrary checksums, assuming that

  • there is a fixed frame length k
  • the last m bits in a frame are a checksum
  • the checksum is formed by a (possibly non-linear) combination of the first k-m bits in the frame
  • there are no collisions between checksums, i.e., that a given sequence of k-m data bits always map onto the same checksum

The test for a T207 checksum looks like this (thanks to Antonio for providing a bit stream):

T207_test.m

Suppose that we only know that the data consist of frames of 14 bits where the last two bits are a checksum but not know how the checksum is formed and where the frames start. Then a search for the presence of such frames results in the plot below using test_checksum.m

test_checksum(bits,14,2)
Indeed the right start position of the frames is recovered when searching for the minimum number of checksum collisions without assuming anything about the way the checksum is formed besides the frame length and the number and position of the checksum bits.

This method can be easily generalized for checksums which are formed using bits both from the present and from the previous frame, e.g., the GPS navigation message, see 20.3.2.5 in IS-GPS-200H.